getwd()
## [1] "/Users/alexg/R files/hair_cortisol/skew-normal FINAL"
library(readxl)
library(psych)
library(dlookr)
## Registered S3 methods overwritten by 'dlookr':
## method from
## plot.transform scales
## print.transform scales
##
## Attaching package: 'dlookr'
## The following object is masked from 'package:psych':
##
## describe
## The following object is masked from 'package:base':
##
## transform
library(vtable)
## Loading required package: kableExtra
library(dplyr)
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:kableExtra':
##
## group_rows
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(reshape)
##
## Attaching package: 'reshape'
## The following object is masked from 'package:dplyr':
##
## rename
library(ggplot2)
##
## Attaching package: 'ggplot2'
## The following objects are masked from 'package:psych':
##
## %+%, alpha
library(brms)
## Loading required package: Rcpp
## Loading 'brms' package (version 2.22.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').
##
## Attaching package: 'brms'
## The following object is masked from 'package:psych':
##
## cs
## The following object is masked from 'package:stats':
##
## ar
library(rethinking)
## Loading required package: cmdstanr
## This is cmdstanr version 0.8.0
## - CmdStanR documentation and vignettes: mc-stan.org/cmdstanr
## - CmdStan path: /Users/alexg/.cmdstan/cmdstan-2.36.0
## - CmdStan version: 2.36.0
## Loading required package: posterior
## This is posterior version 1.6.1
##
## Attaching package: 'posterior'
## The following object is masked from 'package:dlookr':
##
## entropy
## The following objects are masked from 'package:stats':
##
## mad, sd, var
## The following objects are masked from 'package:base':
##
## %in%, match
## Loading required package: parallel
## rethinking (Version 2.42)
##
## Attaching package: 'rethinking'
## The following objects are masked from 'package:brms':
##
## LOO, stancode, WAIC
## The following objects are masked from 'package:psych':
##
## logistic, logit, sim
## The following object is masked from 'package:stats':
##
## rstudent
library(loo)
## This is loo version 2.8.0
## - Online documentation and vignettes at mc-stan.org/loo
## - As of v2.0.0 loo defaults to 1 core but we recommend using as many as possible. Use the 'cores' argument or set options(mc.cores = NUM_CORES) for an entire session.
##
## Attaching package: 'loo'
## The following object is masked from 'package:rethinking':
##
## compare
library(priorsense)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ lubridate 1.9.4 ✔ tibble 3.2.1
## ✔ purrr 1.0.4 ✔ tidyr 1.3.1
## ✔ readr 2.1.5
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ ggplot2::%+%() masks psych::%+%()
## ✖ ggplot2::alpha() masks psych::alpha()
## ✖ tidyr::expand() masks reshape::expand()
## ✖ tidyr::extract() masks dlookr::extract()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::group_rows() masks kableExtra::group_rows()
## ✖ dplyr::lag() masks stats::lag()
## ✖ purrr::map() masks rethinking::map()
## ✖ reshape::rename() masks dplyr::rename()
## ✖ lubridate::stamp() masks reshape::stamp()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(vioplot)
## Loading required package: sm
## Package 'sm', version 2.2-6.0: type help(sm) for summary information
##
## Attaching package: 'sm'
##
## The following object is masked from 'package:dlookr':
##
## binning
##
## Loading required package: zoo
##
## Attaching package: 'zoo'
##
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(bayesplot)
## This is bayesplot version 1.12.0
## - Online documentation and vignettes at mc-stan.org/bayesplot
## - bayesplot theme set to bayesplot::theme_default()
## * Does _not_ affect other ggplot2 plots
## * See ?bayesplot_theme_set for details on theme setting
##
## Attaching package: 'bayesplot'
##
## The following object is masked from 'package:posterior':
##
## rhat
##
## The following object is masked from 'package:brms':
##
## rhat
library(bayestestR)
df <- read_xlsx("hair_cort_dog_all.xlsx", col_types = c("text", "text",
"text", "text", "text", "text",
"text", "numeric","text", "skip",
"numeric", "skip", "skip",
"numeric", "skip"))
df <- as.data.frame(df)
dim(df) # will tell you how many rows and columns the dataset has
## [1] 73 11
class(df) # will tell you what data structure has the dataset been assigned
## [1] "data.frame"
head(df)
## number group visit season breed_group coat_colour sex age comorbidity
## 1 c1 stopped v0 winter ret dark Male 43 yes
## 2 c2 stopped v0 autumn mix dark Male 105 yes
## 3 c3 stopped v0 spring ckcs mix Female 117 yes
## 4 c4 stopped v0 summer ret dark Female 108 yes
## 5 c5 stopped v0 summer ret dark Female 110 yes
## 6 c6 stopped v0 winter mix light Female 120 yes
## fat_percent cortisol
## 1 52.21393 4.924220
## 2 38.52059 7.304202
## 3 46.94916 1.590000
## 4 44.46813 0.861570
## 5 39.59363 6.217317
## 6 NA 4.426785
numeric_df <- Filter(is.numeric, df)
describe(numeric_df) # the describe function in psych provides summary stats
## # A tibble: 3 × 26
## described_variables n na mean sd se_mean IQR skewness kurtosis
## <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 age 73 0 95.8 35.6 4.16 44 -0.104 -0.00589
## 2 fat_percent 55 18 40.5 7.82 1.05 7.82 -0.294 1.12
## 3 cortisol 73 0 8.11 16.5 1.93 5.43 4.05 18.7
## # ℹ 17 more variables: p00 <dbl>, p01 <dbl>, p05 <dbl>, p10 <dbl>, p20 <dbl>,
## # p25 <dbl>, p30 <dbl>, p40 <dbl>, p50 <dbl>, p60 <dbl>, p70 <dbl>,
## # p75 <dbl>, p80 <dbl>, p90 <dbl>, p95 <dbl>, p99 <dbl>, p100 <dbl>
plot_normality(numeric_df)
apply(numeric_df, 2, shapiro.test)
## $age
##
## Shapiro-Wilk normality test
##
## data: newX[, i]
## W = 0.97361, p-value = 0.1288
##
##
## $fat_percent
##
## Shapiro-Wilk normality test
##
## data: newX[, i]
## W = 0.97956, p-value = 0.4692
##
##
## $cortisol
##
## Shapiro-Wilk normality test
##
## data: newX[, i]
## W = 0.46269, p-value = 6.756e-15
qqnorm(df$cortisol)
qqline(df$cortisol, col = "red")
qqnorm(log(df$cortisol))
qqline(log(df$cortisol), col = "red")
shapiro.test(log(df$cortisol))
##
## Shapiro-Wilk normality test
##
## data: log(df$cortisol)
## W = 0.94725, p-value = 0.004126
summary(df$cortisol)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.4141 1.4119 2.3331 8.1089 6.8455 104.6172
summary(log(df$cortisol))
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.8817 0.3449 0.8472 1.1816 1.9236 4.6503
df$lgCort <- log(df$cortisol)
summary(df$lgCort)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.8817 0.3449 0.8472 1.1816 1.9236 4.6503
hist(df$lgCort)
df$breed <- df$breed_group
df$breed <- factor(df$breed, levels = c("mix", "ckcs", "pug", "ret", "other"))
head(df$breed)
## [1] ret mix ckcs ret ret mix
## Levels: mix ckcs pug ret other
df$coat_colour <- factor(df$coat_colour, levels = c("light", "mix", "dark"), ordered = FALSE)
head(df$coat_colour)
## [1] dark dark mix dark dark light
## Levels: light mix dark
sumtable(df)
| Variable | N | Mean | Std. Dev. | Min | Pctl. 25 | Pctl. 75 | Max |
|---|---|---|---|---|---|---|---|
| group | 73 | ||||||
| … completed | 42 | 58% | |||||
| … stopped | 31 | 42% | |||||
| visit | 73 | ||||||
| … v0 | 52 | 71% | |||||
| … v1 | 21 | 29% | |||||
| season | 73 | ||||||
| … autumn | 21 | 29% | |||||
| … spring | 14 | 19% | |||||
| … summer | 22 | 30% | |||||
| … winter | 16 | 22% | |||||
| breed_group | 73 | ||||||
| … ckcs | 7 | 10% | |||||
| … mix | 16 | 22% | |||||
| … other | 26 | 36% | |||||
| … pug | 7 | 10% | |||||
| … ret | 17 | 23% | |||||
| coat_colour | 73 | ||||||
| … light | 27 | 37% | |||||
| … mix | 16 | 22% | |||||
| … dark | 30 | 41% | |||||
| sex | 73 | ||||||
| … Female | 43 | 59% | |||||
| … Male | 30 | 41% | |||||
| age | 73 | 96 | 36 | 16 | 73 | 117 | 182 |
| comorbidity | 73 | ||||||
| … no | 15 | 21% | |||||
| … yes | 58 | 79% | |||||
| fat_percent | 55 | 40 | 7.8 | 18 | 37 | 45 | 61 |
| cortisol | 73 | 8.1 | 16 | 0.41 | 1.4 | 6.8 | 105 |
| lgCort | 73 | 1.2 | 1.2 | -0.88 | 0.34 | 1.9 | 4.7 |
| breed | 73 | ||||||
| … mix | 16 | 22% | |||||
| … ckcs | 7 | 10% | |||||
| … pug | 7 | 10% | |||||
| … ret | 17 | 23% | |||||
| … other | 26 | 36% |
df_stopped <- df[df$group == "stopped", ]
sumtable(df_stopped, add.median = TRUE)
| Variable | N | Mean | Std. Dev. | Min | Pctl. 25 | Pctl. 50 | Pctl. 75 | Max |
|---|---|---|---|---|---|---|---|---|
| group | 31 | |||||||
| … stopped | 31 | 100% | ||||||
| visit | 31 | |||||||
| … v0 | 31 | 100% | ||||||
| season | 31 | |||||||
| … autumn | 10 | 32% | ||||||
| … spring | 9 | 29% | ||||||
| … summer | 5 | 16% | ||||||
| … winter | 7 | 23% | ||||||
| breed_group | 31 | |||||||
| … ckcs | 5 | 16% | ||||||
| … mix | 6 | 19% | ||||||
| … other | 10 | 32% | ||||||
| … pug | 3 | 10% | ||||||
| … ret | 7 | 23% | ||||||
| coat_colour | 31 | |||||||
| … light | 8 | 26% | ||||||
| … mix | 10 | 32% | ||||||
| … dark | 13 | 42% | ||||||
| sex | 31 | |||||||
| … Female | 17 | 55% | ||||||
| … Male | 14 | 45% | ||||||
| age | 31 | 90 | 34 | 16 | 71 | 98 | 110 | 155 |
| comorbidity | 31 | |||||||
| … no | 1 | 3% | ||||||
| … yes | 30 | 97% | ||||||
| fat_percent | 26 | 41 | 5.5 | 31 | 38 | 42 | 44 | 54 |
| cortisol | 31 | 11 | 22 | 0.46 | 1.8 | 2.8 | 7.2 | 105 |
| lgCort | 31 | 1.4 | 1.3 | -0.77 | 0.57 | 1 | 2 | 4.7 |
| breed | 31 | |||||||
| … mix | 6 | 19% | ||||||
| … ckcs | 5 | 16% | ||||||
| … pug | 3 | 10% | ||||||
| … ret | 7 | 23% | ||||||
| … other | 10 | 32% |
df_completed <- df[df$group == "completed", ]
df_compv0 <- df_completed[df_completed$visit == "v0", ]
df_compv1 <- df_completed[df_completed$visit == "v1", ]
sumtable(df_compv0, add.median = TRUE)
| Variable | N | Mean | Std. Dev. | Min | Pctl. 25 | Pctl. 50 | Pctl. 75 | Max |
|---|---|---|---|---|---|---|---|---|
| group | 21 | |||||||
| … completed | 21 | 100% | ||||||
| visit | 21 | |||||||
| … v0 | 21 | 100% | ||||||
| season | 21 | |||||||
| … autumn | 3 | 14% | ||||||
| … spring | 4 | 19% | ||||||
| … summer | 8 | 38% | ||||||
| … winter | 6 | 29% | ||||||
| breed_group | 21 | |||||||
| … ckcs | 1 | 5% | ||||||
| … mix | 5 | 24% | ||||||
| … other | 8 | 38% | ||||||
| … pug | 2 | 10% | ||||||
| … ret | 5 | 24% | ||||||
| coat_colour | 21 | |||||||
| … light | 10 | 48% | ||||||
| … mix | 3 | 14% | ||||||
| … dark | 8 | 38% | ||||||
| sex | 21 | |||||||
| … Female | 13 | 62% | ||||||
| … Male | 8 | 38% | ||||||
| age | 21 | 91 | 36 | 29 | 67 | 96 | 109 | 165 |
| comorbidity | 21 | |||||||
| … no | 7 | 33% | ||||||
| … yes | 14 | 67% | ||||||
| fat_percent | 17 | 45 | 6.4 | 37 | 40 | 43 | 48 | 61 |
| cortisol | 21 | 7.6 | 14 | 0.69 | 1.7 | 2.2 | 6.8 | 59 |
| lgCort | 21 | 1.1 | 1.2 | -0.38 | 0.52 | 0.77 | 1.9 | 4.1 |
| breed | 21 | |||||||
| … mix | 5 | 24% | ||||||
| … ckcs | 1 | 5% | ||||||
| … pug | 2 | 10% | ||||||
| … ret | 5 | 24% | ||||||
| … other | 8 | 38% |
sumtable(df_compv1, add.median = TRUE)
| Variable | N | Mean | Std. Dev. | Min | Pctl. 25 | Pctl. 50 | Pctl. 75 | Max |
|---|---|---|---|---|---|---|---|---|
| group | 21 | |||||||
| … completed | 21 | 100% | ||||||
| visit | 21 | |||||||
| … v1 | 21 | 100% | ||||||
| season | 21 | |||||||
| … autumn | 8 | 38% | ||||||
| … spring | 1 | 5% | ||||||
| … summer | 9 | 43% | ||||||
| … winter | 3 | 14% | ||||||
| breed_group | 21 | |||||||
| … ckcs | 1 | 5% | ||||||
| … mix | 5 | 24% | ||||||
| … other | 8 | 38% | ||||||
| … pug | 2 | 10% | ||||||
| … ret | 5 | 24% | ||||||
| coat_colour | 21 | |||||||
| … light | 9 | 43% | ||||||
| … mix | 3 | 14% | ||||||
| … dark | 9 | 43% | ||||||
| sex | 21 | |||||||
| … Female | 13 | 62% | ||||||
| … Male | 8 | 38% | ||||||
| age | 21 | 109 | 35 | 42 | 94 | 105 | 129 | 182 |
| comorbidity | 21 | |||||||
| … no | 7 | 33% | ||||||
| … yes | 14 | 67% | ||||||
| fat_percent | 12 | 33 | 8.8 | 18 | 27 | 33 | 38 | 50 |
| cortisol | 21 | 3.6 | 3.7 | 0.41 | 1.2 | 1.8 | 5.6 | 15 |
| lgCort | 21 | 0.85 | 0.94 | -0.88 | 0.21 | 0.58 | 1.7 | 2.7 |
| breed | 21 | |||||||
| … mix | 5 | 24% | ||||||
| … ckcs | 1 | 5% | ||||||
| … pug | 2 | 10% | ||||||
| … ret | 5 | 24% | ||||||
| … other | 8 | 38% |
## three way cross tabs (xtabs) and flatten the table
ftable(xtabs(~ visit + group, data = df))
## group completed stopped
## visit
## v0 21 31
## v1 21 0
par(mfrow = c(1,1))
vioplot(cortisol ~ sex, col = "firebrick",
data = df)
par(mfrow = c(1,1))
vioplot(lgCort ~ sex, col = "lemonchiffon",
data = df)
par(mfrow = c(1,1))
vioplot(lgCort ~ breed, col = "firebrick",
data = df)
stripchart(lgCort ~ breed, vertical = TRUE, method = "jitter",
col = "steelblue3", data = df, pch = 20)
par(mfrow = c(1,1))
vioplot(cortisol ~ comorbidity, col = "steelblue",
data = df)
par(mfrow = c(1,1))
vioplot(lgCort ~ comorbidity, col = "steelblue",
data = df)
plot(cortisol ~ age, col = "steelblue3", data = df, pch = 20)
plot(lgCort ~ age, col = "steelblue3", data = df, pch = 20)
plot(cortisol ~ fat_percent, col = "steelblue3", data = df, pch = 20)
plot(lgCort ~ fat_percent, col = "steelblue3", data = df, pch = 20)
corr.test(df$lgCort, df$fat_percent)
## Call:corr.test(x = df$lgCort, y = df$fat_percent)
## Correlation matrix
## [1] 0.06
## Sample Size
## [1] 55
## These are the unadjusted probability values.
## The probability values adjusted for multiple tests are in the p.adj object.
## [1] 0.65
##
## To see confidence intervals of the correlations, print with the short=FALSE option
df$sFat <- standardize(df$fat_percent)
summary(df$sFat)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## -2.89165 -0.43568 0.04814 0.00000 0.56366 2.64873 18
hist(df$sFat)
df$sAge <- standardize(df$age)
summary(df$sAge)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2.2444 -0.6422 0.1448 0.0000 0.5945 2.4215
df$slgCort <- standardize(df$lgC)
summary(df$slgCort)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.7079 -0.6925 -0.2768 0.0000 0.6142 2.8713
hist(df$slgCort, breaks =20, col = "steelblue3", main = "(a) Histogram of observed log hair cortisol", xlab = "Log hair cortisol (standardised)", xlim = c(-2, 3))
x <- seq(-2, 4, length.out = 100)
y <- dskew_normal(x, mu = 0, sigma = 1, alpha = 4)
plot(y ~ x, type = "l", main = "(b) Skew-normal distribution with alpha of 4", ylab = "Density", xlab = "Log hair cortisol (standardised)", xlim = c(-2, 3), lwd = 2, col = "steelblue3")
par(mfrow = c(1,2))
hist(df$slgCort, breaks =20, col = "steelblue3", main = "(a) Log hair cortisol (observed)", xlab = "Log hair cortisol (standardised)", xlim = c(-2, 3))
x <- seq(-2, 4, length.out = 100)
y <- dskew_normal(x, mu = 0, sigma = 1, alpha = 4)
plot(y ~ x, type = "l", main = "(b) Skew-normal distribution (alpha=4)", ylab = "Density", xlab = "Log hair cortisol (standardised)", xlim = c(-2, 3), lwd = 2, col = "steelblue3")
par(mfrow = c(1,1))
# prepare histogram data
hist_data <- df$slgCort
# Prepare skew-normal curve dats
x <- seq(-3, 3, length.out = 100)
y <- dskew_normal(x, mu = 0, sigma = 1, alpha = 4)
# Histogram
hist(hist_data, prob = TRUE, col = "steelblue2", breaks = 20,
ylim = c(0, 0.8), # Adjust y-axis limits to fit the data
xlim = c(-3, 3),
main = "Observed hair cortisol vs skew-normal",
xlab = "Log hair cortisol (standardised)",
ylab = "Density",
font.lab = 2)
# Add skew-normal curve
lines(x, y, col = "firebrick3", lwd = 3)
# Add normal distribution for comparison
lines(x, dnorm(x, mean = 0, sd = 1), col = "darkblue", lwd = 2, lty = 2)
x <- seq(-3, 3, length.out = 100)
y <- dnorm(x, mean = 0, sd = 0.5)
plot(y ~ x, type = "l",
col = "steelblue3", lwd = 3,
main = "Prior for the effect of body fat on log hair cortisol",
xlab = "Prior probability distribution",
ylab = "Density",
font.lab = 2)
abline(v = 0, lwd = 0.5)
abline(v = 0, col = "firebrick3", lty = 2, lwd = 2)
df2 <- na.omit(df)
model <- brm(slgCort ~ sFat + sAge + breed + sex + comorbidity + (1 | visit), family = skew_normal(), data = df2)
default_prior(slgCort ~ sFat + sAge + breed + sex + comorbidity + (1 | visit),
family = skew_normal(),
data = df2)
## prior class coef group resp dpar nlpar lb ub
## normal(0, 4) alpha
## (flat) b
## (flat) b breedckcs
## (flat) b breedother
## (flat) b breedpug
## (flat) b breedret
## (flat) b comorbidityyes
## (flat) b sAge
## (flat) b sexMale
## (flat) b sFat
## student_t(3, -0.1, 2.5) Intercept
## student_t(3, 0, 2.5) sd 0
## student_t(3, 0, 2.5) sd visit 0
## student_t(3, 0, 2.5) sd Intercept visit 0
## student_t(3, 0, 2.5) sigma 0
## source
## default
## default
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## default
## default
## (vectorized)
## (vectorized)
## default
exp(0.983) # obese
## [1] 2.672462
exp(0.862) # normal
## [1] 2.367892
exp(0.904) # overweight
## [1] 2.469461
Change from normal to obese = 2.37 to 2.67, so only ~10% increase across the range. Prior allowing beta increase of 0.1 per SD increase of BF would work,meaning that we expect change to be between -0.1 and +0.3 Include regularising function to improve efficiency
Change in cortisol when slgCort increases from mean (0) by 0.1 X 1 SD ~10% increase in cortisol for ~10% increase in fat
mlgC <- mean(df2$lgCort)
sdlgC <- sd(df2$lgCort)
exp(mlgC)
## [1] 3.600573
exp(mlgC + 0.1 * sdlgC)
## [1] 4.083976
Change in fat when increases from mean by 1 SD 10% increase in fat
mF <- mean(df2$fat_percent)
sdF <- sd(df2$fat_percent)
mF
## [1] 40.47973
mF + sdF
## [1] 48.30141
Suggests that a 10% increase in fat would be associated with a 10% increase in mean cortisol
Against this evidence in dogs, BCS had negative impact on log hair cortisol… albeit for dogs with BCS 1-6.
NB Bowland found no age effect. They also found a negative effect of BCS on log hair cortisol (beta -0.03). However, BCS ranged from 1-6 (with only 1 BCS 6), so few overweight…. could suggest poor nutrition and health cf obesity. (Bowland JB. Front. Vet. Sci 2020; 7:565346. doi: 10.3389/fvets.2020.565346)
Nonetheless, safer to use a neutral regularising prior.
hair_cort is greater in dogs with atopy vs control Ref: Vet
Dermatol 2023 Aug;34(4):339-347. doi: 10.1111/vde.13151.
hair cortisol before treatment: 0.605 hair cortisol after treatment:
0.25 treatment leads to ~50% reduction
hair_cort greater when dogs are stressed (shelter dogs) Ref: Scientific Reports | (2022) 12:5117 | https://doi.org/10.1038/s41598-022-09140-w t0 16.0 +/- 6.8 t1 21.8 +/- 9.4 Average increase = 30-40%
hair_cort greater when behavioural stress (competition) and seasonal affect (but variable) Ref: Scientific Reports | 6:19631 | DOI: 10.1038/srep19631 Greater cortisol in January, although only in competing dogs. Cortisol > in competing cf companion and working 40 pg/mg vs 10-20, so approx double.
hair_cort greater when behavioural stress (agility) Ref: PLOS ONE | https://doi.org/10.1371/journal.pone.0216000 yes 0.26 (0.17-0.35) not 0.15 (0.10-0.27) Average increase = 50%
Overall, as not exactly the same as effect of various diseases, set a modest (and regularising) prior… 20% increase in cortisol (0.2 SD) for presence of comorbidity
In humans, males have > hair cortisol cf females (Binz TM. ForensicSciInt 2018; 284:33–8. doi: 10.1016/j.forsciint.2017.12.032) …but effect is opposite in vervet monlkys (Laudenslager ML. Psychoneuroendocrinology 2012; 37:1736–9. doi: 10.1016/j.psyneuen.2012.03.015). … no effect of sex in a previous dog study (Macbeth BJ. Wildl Soc Bull 2012; 36:747–58. doi: 10.1002/wsb.219) … however, this study did find that neutered dogs had decreased hair cortisol. … as all dogs in the study were neutered, this means that an effect of sex is less likely. Bowland et al, female dogs had > hair cortisol than male dogs, but all dogs were intact (Bowland JB. Front. Vet. Sci 2020; 7:565346. doi: 10.3389/fvets.2020.565346)
Further, this effect was lost when accounting for other effects.
Therefore, use a regularising prior but keep it neutral and broad, to allow the effect to be either way.
NB Bowland found no age effect. They also found a negative effect of BCS on log hair cortisol (beta -0.03). However, BCS ranged from 1-6 (with only 1 BCS 6), so few overweight…. could suggest poor nutrition and health cf obesity. (Bowland JB. Front. Vet. Sci 2020; 7:565346. doi: 10.3389/fvets.2020.565346)
No published data about effects of breed on hair cortisol, but this is plausible However, unclear as to which breeds will differ and which way. Therefore, use a regularising prior but keep it neutral and broad.
# Set individual priors
prior_int <- set_prior("normal(0, 0.5)", class = "Intercept")
prior_b <- set_prior("normal(0, 1)", class = "b") # betas for breed
prior_b_sex <- set_prior("normal(0, 1)", class = "b", coef = "sexMale")
prior_b_co <- set_prior("normal(0.25, 1)", class = "b", coef = "comorbidityyes")
prior_b_f <- set_prior("normal(0, 0.5)", class = "b", coef = "sFat")
prior_b_sAge <- set_prior("normal(0, 0.5)", class = "b", coef = "sAge")
prior_sd <- set_prior("normal(0, 1)", class = "sd")
prior_alpha <- set_prior("normal(4, 2)", class = "alpha")
# Combine priors into list
priors <- c(prior_int, prior_b, prior_b_sex, prior_b_co, prior_b_f, prior_b_sAge, prior_sd, prior_alpha)
NB, as we have comorbidity in the nodek, we allowing sigma to vary across comorbidity.
par(mfrow = c(1,1))
x <- seq(-3, 3, length.out = 100)
y <- dnorm(x, mean = 0, sd = 0.5)
plot(y ~ x, type = "l")
x <- seq(0, 3, length.out = 100)
y <- dexp(x, rate = 1)
plot(y ~ x, type = "l")
Based on distribution of log normal hair cortisol, expect things to be skewed to the right. Try different levels of alpha for skew normal.
x <- seq(-2, 4, length.out = 100)
y <- dskew_normal(x, mu = 0, sigma = 1, alpha = 4)
plot(y ~ x, type = "l")
Given we have a right skew, alpha needs to be positive and 4 appears to reflect the skew seen in the data… keep sd briad (eg 2) broad to allow some flexibility
x <- seq(0, 5, length.out = 100)
y <- dnorm(x, mean = 0, sd = 1)
plot(y ~ x, type = "l")
x <- seq(-3, 3, length.out = 100)
y <- dnorm(x, mean = 0, sd = 1.0)
plot(y ~ x, type = "l")
x <- seq(-3, 3, length.out = 100)
y <- dnorm(x, mean = 0, sd = 0.5)
plot(y ~ x, type = "l")
x <- seq(-3, 3, length.out = 100)
y <- dnorm(x, mean = 0.25, sd =1)
plot(y ~ x, type = "l", col = "steelblue3", lwd = 3,
main = "Prior for comorbidity effect on log hair cortisol",
xlab = "Prior probability distribution",
ylab = "Density",
font.lab = 2)
abline(v = 0, lwd = 0.5)
abline(v = 0.25, col = "firebrick3", lty = 2, lwd = 2)
x <- seq(-3, 3, length.out = 100)
y <- dnorm(x, mean = 0, sd = 1)
plot(y ~ x, type = "l")
Increased adapt_delta >0.8 (0.9999 here), as had divergent transitions
set.seed(666)
model <- brm(bf(slgCort ~ sFat + sAge + breed + sex + comorbidity + (1 | visit),
sigma ~ comorbidity),
family = skew_normal(),
prior = priors,
data = df2,
control=list(adapt_delta=0.999999, stepsize = 0.001, max_treedepth =15),
iter = 8000, warmup = 2000,
cores = 4,
save_pars = save_pars(all =TRUE),
sample_prior = TRUE)
## Compiling Stan program...
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 17.0.0 (clang-1700.0.13.5)’
## using SDK: ‘MacOSX15.5.sdk’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/opt/R/arm64/include -fPIC -falign-functions=64 -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
## 679 | #include <cmath>
## | ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1
## Start sampling
## Found more than one class "stanfit" in cache; using the first, from namespace 'rethinking'
## Also defined by 'rstan'
## Found more than one class "stanfit" in cache; using the first, from namespace 'rethinking'
## Also defined by 'rstan'
## Found more than one class "stanfit" in cache; using the first, from namespace 'rethinking'
## Also defined by 'rstan'
## Found more than one class "stanfit" in cache; using the first, from namespace 'rethinking'
## Also defined by 'rstan'
## Warning: There were 3 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: Examine the pairs() plot to diagnose sampling problems
summary(model)
## Warning: There were 3 divergent transitions after warmup. Increasing
## adapt_delta above 0.999999 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Family: skew_normal
## Links: mu = identity; sigma = log; alpha = identity
## Formula: slgCort ~ sFat + sAge + breed + sex + comorbidity + (1 | visit)
## sigma ~ comorbidity
## Data: df2 (Number of observations: 55)
## Draws: 4 chains, each with iter = 8000; warmup = 2000; thin = 1;
## total post-warmup draws = 24000
##
## Multilevel Hyperparameters:
## ~visit (Number of levels: 2)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.39 0.37 0.01 1.40 1.00 10046 10684
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept -0.57 0.43 -1.44 0.26 1.00 11140
## sigma_Intercept -0.94 0.34 -1.52 -0.19 1.00 11614
## sFat 0.33 0.13 0.04 0.56 1.00 14598
## sAge 0.06 0.13 -0.21 0.32 1.00 20741
## breedckcs -0.11 0.44 -0.99 0.73 1.00 19374
## breedpug -0.17 0.44 -1.03 0.71 1.00 13867
## breedret -0.23 0.33 -0.88 0.43 1.00 14992
## breedother 0.13 0.34 -0.52 0.83 1.00 13159
## sexMale 0.10 0.24 -0.38 0.58 1.00 20484
## comorbidityyes 0.80 0.25 0.29 1.28 1.00 16599
## sigma_comorbidityyes 1.09 0.36 0.31 1.73 1.00 12298
## Tail_ESS
## Intercept 12627
## sigma_Intercept 12112
## sFat 14824
## sAge 17258
## breedckcs 17407
## breedpug 13352
## breedret 14849
## breedother 14221
## sexMale 16305
## comorbidityyes 14612
## sigma_comorbidityyes 12125
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## alpha 4.22 1.66 1.19 7.69 1.00 18797 16815
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
color_scheme_set("blue")
plot(model)
Looking for hairy caterpillars
mcmc_plot(model, type = 'rank_overlay')
posterior <- as.array(model)
dimnames(posterior)
## $iteration
## [1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10"
## [11] "11" "12" "13" "14" "15" "16" "17" "18" "19" "20"
## [21] "21" "22" "23" "24" "25" "26" "27" "28" "29" "30"
## [31] "31" "32" "33" "34" "35" "36" "37" "38" "39" "40"
## [41] "41" "42" "43" "44" "45" "46" "47" "48" "49" "50"
## [51] "51" "52" "53" "54" "55" "56" "57" "58" "59" "60"
## [61] "61" "62" "63" "64" "65" "66" "67" "68" "69" "70"
## [71] "71" "72" "73" "74" "75" "76" "77" "78" "79" "80"
## [81] "81" "82" "83" "84" "85" "86" "87" "88" "89" "90"
## [91] "91" "92" "93" "94" "95" "96" "97" "98" "99" "100"
## [101] "101" "102" "103" "104" "105" "106" "107" "108" "109" "110"
## [111] "111" "112" "113" "114" "115" "116" "117" "118" "119" "120"
## [121] "121" "122" "123" "124" "125" "126" "127" "128" "129" "130"
## [131] "131" "132" "133" "134" "135" "136" "137" "138" "139" "140"
## [141] "141" "142" "143" "144" "145" "146" "147" "148" "149" "150"
## [151] "151" "152" "153" "154" "155" "156" "157" "158" "159" "160"
## [161] "161" "162" "163" "164" "165" "166" "167" "168" "169" "170"
## [171] "171" "172" "173" "174" "175" "176" "177" "178" "179" "180"
## [181] "181" "182" "183" "184" "185" "186" "187" "188" "189" "190"
## [191] "191" "192" "193" "194" "195" "196" "197" "198" "199" "200"
## [201] "201" "202" "203" "204" "205" "206" "207" "208" "209" "210"
## [211] "211" "212" "213" "214" "215" "216" "217" "218" "219" "220"
## [221] "221" "222" "223" "224" "225" "226" "227" "228" "229" "230"
## [231] "231" "232" "233" "234" "235" "236" "237" "238" "239" "240"
## [241] "241" "242" "243" "244" "245" "246" "247" "248" "249" "250"
## [251] "251" "252" "253" "254" "255" "256" "257" "258" "259" "260"
## [261] "261" "262" "263" "264" "265" "266" "267" "268" "269" "270"
## [271] "271" "272" "273" "274" "275" "276" "277" "278" "279" "280"
## [281] "281" "282" "283" "284" "285" "286" "287" "288" "289" "290"
## [291] "291" "292" "293" "294" "295" "296" "297" "298" "299" "300"
## [301] "301" "302" "303" "304" "305" "306" "307" "308" "309" "310"
## [311] "311" "312" "313" "314" "315" "316" "317" "318" "319" "320"
## [321] "321" "322" "323" "324" "325" "326" "327" "328" "329" "330"
## [331] "331" "332" "333" "334" "335" "336" "337" "338" "339" "340"
## [341] "341" "342" "343" "344" "345" "346" "347" "348" "349" "350"
## [351] "351" "352" "353" "354" "355" "356" "357" "358" "359" "360"
## [361] "361" "362" "363" "364" "365" "366" "367" "368" "369" "370"
## [371] "371" "372" "373" "374" "375" "376" "377" "378" "379" "380"
## [381] "381" "382" "383" "384" "385" "386" "387" "388" "389" "390"
## [391] "391" "392" "393" "394" "395" "396" "397" "398" "399" "400"
## [401] "401" "402" "403" "404" "405" "406" "407" "408" "409" "410"
## [411] "411" "412" "413" "414" "415" "416" "417" "418" "419" "420"
## [421] "421" "422" "423" "424" "425" "426" "427" "428" "429" "430"
## [431] "431" "432" "433" "434" "435" "436" "437" "438" "439" "440"
## [441] "441" "442" "443" "444" "445" "446" "447" "448" "449" "450"
## [451] "451" "452" "453" "454" "455" "456" "457" "458" "459" "460"
## [461] "461" "462" "463" "464" "465" "466" "467" "468" "469" "470"
## [471] "471" "472" "473" "474" "475" "476" "477" "478" "479" "480"
## [481] "481" "482" "483" "484" "485" "486" "487" "488" "489" "490"
## [491] "491" "492" "493" "494" "495" "496" "497" "498" "499" "500"
## [501] "501" "502" "503" "504" "505" "506" "507" "508" "509" "510"
## [511] "511" "512" "513" "514" "515" "516" "517" "518" "519" "520"
## [521] "521" "522" "523" "524" "525" "526" "527" "528" "529" "530"
## [531] "531" "532" "533" "534" "535" "536" "537" "538" "539" "540"
## [541] "541" "542" "543" "544" "545" "546" "547" "548" "549" "550"
## [551] "551" "552" "553" "554" "555" "556" "557" "558" "559" "560"
## [561] "561" "562" "563" "564" "565" "566" "567" "568" "569" "570"
## [571] "571" "572" "573" "574" "575" "576" "577" "578" "579" "580"
## [581] "581" "582" "583" "584" "585" "586" "587" "588" "589" "590"
## [591] "591" "592" "593" "594" "595" "596" "597" "598" "599" "600"
## [601] "601" "602" "603" "604" "605" "606" "607" "608" "609" "610"
## [611] "611" "612" "613" "614" "615" "616" "617" "618" "619" "620"
## [621] "621" "622" "623" "624" "625" "626" "627" "628" "629" "630"
## [631] "631" "632" "633" "634" "635" "636" "637" "638" "639" "640"
## [641] "641" "642" "643" "644" "645" "646" "647" "648" "649" "650"
## [651] "651" "652" "653" "654" "655" "656" "657" "658" "659" "660"
## [661] "661" "662" "663" "664" "665" "666" "667" "668" "669" "670"
## [671] "671" "672" "673" "674" "675" "676" "677" "678" "679" "680"
## [681] "681" "682" "683" "684" "685" "686" "687" "688" "689" "690"
## [691] "691" "692" "693" "694" "695" "696" "697" "698" "699" "700"
## [701] "701" "702" "703" "704" "705" "706" "707" "708" "709" "710"
## [711] "711" "712" "713" "714" "715" "716" "717" "718" "719" "720"
## [721] "721" "722" "723" "724" "725" "726" "727" "728" "729" "730"
## [731] "731" "732" "733" "734" "735" "736" "737" "738" "739" "740"
## [741] "741" "742" "743" "744" "745" "746" "747" "748" "749" "750"
## [751] "751" "752" "753" "754" "755" "756" "757" "758" "759" "760"
## [761] "761" "762" "763" "764" "765" "766" "767" "768" "769" "770"
## [771] "771" "772" "773" "774" "775" "776" "777" "778" "779" "780"
## [781] "781" "782" "783" "784" "785" "786" "787" "788" "789" "790"
## [791] "791" "792" "793" "794" "795" "796" "797" "798" "799" "800"
## [801] "801" "802" "803" "804" "805" "806" "807" "808" "809" "810"
## [811] "811" "812" "813" "814" "815" "816" "817" "818" "819" "820"
## [821] "821" "822" "823" "824" "825" "826" "827" "828" "829" "830"
## [831] "831" "832" "833" "834" "835" "836" "837" "838" "839" "840"
## [841] "841" "842" "843" "844" "845" "846" "847" "848" "849" "850"
## [851] "851" "852" "853" "854" "855" "856" "857" "858" "859" "860"
## [861] "861" "862" "863" "864" "865" "866" "867" "868" "869" "870"
## [871] "871" "872" "873" "874" "875" "876" "877" "878" "879" "880"
## [881] "881" "882" "883" "884" "885" "886" "887" "888" "889" "890"
## [891] "891" "892" "893" "894" "895" "896" "897" "898" "899" "900"
## [901] "901" "902" "903" "904" "905" "906" "907" "908" "909" "910"
## [911] "911" "912" "913" "914" "915" "916" "917" "918" "919" "920"
## [921] "921" "922" "923" "924" "925" "926" "927" "928" "929" "930"
## [931] "931" "932" "933" "934" "935" "936" "937" "938" "939" "940"
## [941] "941" "942" "943" "944" "945" "946" "947" "948" "949" "950"
## [951] "951" "952" "953" "954" "955" "956" "957" "958" "959" "960"
## [961] "961" "962" "963" "964" "965" "966" "967" "968" "969" "970"
## [971] "971" "972" "973" "974" "975" "976" "977" "978" "979" "980"
## [981] "981" "982" "983" "984" "985" "986" "987" "988" "989" "990"
## [991] "991" "992" "993" "994" "995" "996" "997" "998" "999" "1000"
## [1001] "1001" "1002" "1003" "1004" "1005" "1006" "1007" "1008" "1009" "1010"
## [1011] "1011" "1012" "1013" "1014" "1015" "1016" "1017" "1018" "1019" "1020"
## [1021] "1021" "1022" "1023" "1024" "1025" "1026" "1027" "1028" "1029" "1030"
## [1031] "1031" "1032" "1033" "1034" "1035" "1036" "1037" "1038" "1039" "1040"
## [1041] "1041" "1042" "1043" "1044" "1045" "1046" "1047" "1048" "1049" "1050"
## [1051] "1051" "1052" "1053" "1054" "1055" "1056" "1057" "1058" "1059" "1060"
## [1061] "1061" "1062" "1063" "1064" "1065" "1066" "1067" "1068" "1069" "1070"
## [1071] "1071" "1072" "1073" "1074" "1075" "1076" "1077" "1078" "1079" "1080"
## [1081] "1081" "1082" "1083" "1084" "1085" "1086" "1087" "1088" "1089" "1090"
## [1091] "1091" "1092" "1093" "1094" "1095" "1096" "1097" "1098" "1099" "1100"
## [1101] "1101" "1102" "1103" "1104" "1105" "1106" "1107" "1108" "1109" "1110"
## [1111] "1111" "1112" "1113" "1114" "1115" "1116" "1117" "1118" "1119" "1120"
## [1121] "1121" "1122" "1123" "1124" "1125" "1126" "1127" "1128" "1129" "1130"
## [1131] "1131" "1132" "1133" "1134" "1135" "1136" "1137" "1138" "1139" "1140"
## [1141] "1141" "1142" "1143" "1144" "1145" "1146" "1147" "1148" "1149" "1150"
## [1151] "1151" "1152" "1153" "1154" "1155" "1156" "1157" "1158" "1159" "1160"
## [1161] "1161" "1162" "1163" "1164" "1165" "1166" "1167" "1168" "1169" "1170"
## [1171] "1171" "1172" "1173" "1174" "1175" "1176" "1177" "1178" "1179" "1180"
## [1181] "1181" "1182" "1183" "1184" "1185" "1186" "1187" "1188" "1189" "1190"
## [1191] "1191" "1192" "1193" "1194" "1195" "1196" "1197" "1198" "1199" "1200"
## [1201] "1201" "1202" "1203" "1204" "1205" "1206" "1207" "1208" "1209" "1210"
## [1211] "1211" "1212" "1213" "1214" "1215" "1216" "1217" "1218" "1219" "1220"
## [1221] "1221" "1222" "1223" "1224" "1225" "1226" "1227" "1228" "1229" "1230"
## [1231] "1231" "1232" "1233" "1234" "1235" "1236" "1237" "1238" "1239" "1240"
## [1241] "1241" "1242" "1243" "1244" "1245" "1246" "1247" "1248" "1249" "1250"
## [1251] "1251" "1252" "1253" "1254" "1255" "1256" "1257" "1258" "1259" "1260"
## [1261] "1261" "1262" "1263" "1264" "1265" "1266" "1267" "1268" "1269" "1270"
## [1271] "1271" "1272" "1273" "1274" "1275" "1276" "1277" "1278" "1279" "1280"
## [1281] "1281" "1282" "1283" "1284" "1285" "1286" "1287" "1288" "1289" "1290"
## [1291] "1291" "1292" "1293" "1294" "1295" "1296" "1297" "1298" "1299" "1300"
## [1301] "1301" "1302" "1303" "1304" "1305" "1306" "1307" "1308" "1309" "1310"
## [1311] "1311" "1312" "1313" "1314" "1315" "1316" "1317" "1318" "1319" "1320"
## [1321] "1321" "1322" "1323" "1324" "1325" "1326" "1327" "1328" "1329" "1330"
## [1331] "1331" "1332" "1333" "1334" "1335" "1336" "1337" "1338" "1339" "1340"
## [1341] "1341" "1342" "1343" "1344" "1345" "1346" "1347" "1348" "1349" "1350"
## [1351] "1351" "1352" "1353" "1354" "1355" "1356" "1357" "1358" "1359" "1360"
## [1361] "1361" "1362" "1363" "1364" "1365" "1366" "1367" "1368" "1369" "1370"
## [1371] "1371" "1372" "1373" "1374" "1375" "1376" "1377" "1378" "1379" "1380"
## [1381] "1381" "1382" "1383" "1384" "1385" "1386" "1387" "1388" "1389" "1390"
## [1391] "1391" "1392" "1393" "1394" "1395" "1396" "1397" "1398" "1399" "1400"
## [1401] "1401" "1402" "1403" "1404" "1405" "1406" "1407" "1408" "1409" "1410"
## [1411] "1411" "1412" "1413" "1414" "1415" "1416" "1417" "1418" "1419" "1420"
## [1421] "1421" "1422" "1423" "1424" "1425" "1426" "1427" "1428" "1429" "1430"
## [1431] "1431" "1432" "1433" "1434" "1435" "1436" "1437" "1438" "1439" "1440"
## [1441] "1441" "1442" "1443" "1444" "1445" "1446" "1447" "1448" "1449" "1450"
## [1451] "1451" "1452" "1453" "1454" "1455" "1456" "1457" "1458" "1459" "1460"
## [1461] "1461" "1462" "1463" "1464" "1465" "1466" "1467" "1468" "1469" "1470"
## [1471] "1471" "1472" "1473" "1474" "1475" "1476" "1477" "1478" "1479" "1480"
## [1481] "1481" "1482" "1483" "1484" "1485" "1486" "1487" "1488" "1489" "1490"
## [1491] "1491" "1492" "1493" "1494" "1495" "1496" "1497" "1498" "1499" "1500"
## [1501] "1501" "1502" "1503" "1504" "1505" "1506" "1507" "1508" "1509" "1510"
## [1511] "1511" "1512" "1513" "1514" "1515" "1516" "1517" "1518" "1519" "1520"
## [1521] "1521" "1522" "1523" "1524" "1525" "1526" "1527" "1528" "1529" "1530"
## [1531] "1531" "1532" "1533" "1534" "1535" "1536" "1537" "1538" "1539" "1540"
## [1541] "1541" "1542" "1543" "1544" "1545" "1546" "1547" "1548" "1549" "1550"
## [1551] "1551" "1552" "1553" "1554" "1555" "1556" "1557" "1558" "1559" "1560"
## [1561] "1561" "1562" "1563" "1564" "1565" "1566" "1567" "1568" "1569" "1570"
## [1571] "1571" "1572" "1573" "1574" "1575" "1576" "1577" "1578" "1579" "1580"
## [1581] "1581" "1582" "1583" "1584" "1585" "1586" "1587" "1588" "1589" "1590"
## [1591] "1591" "1592" "1593" "1594" "1595" "1596" "1597" "1598" "1599" "1600"
## [1601] "1601" "1602" "1603" "1604" "1605" "1606" "1607" "1608" "1609" "1610"
## [1611] "1611" "1612" "1613" "1614" "1615" "1616" "1617" "1618" "1619" "1620"
## [1621] "1621" "1622" "1623" "1624" "1625" "1626" "1627" "1628" "1629" "1630"
## [1631] "1631" "1632" "1633" "1634" "1635" "1636" "1637" "1638" "1639" "1640"
## [1641] "1641" "1642" "1643" "1644" "1645" "1646" "1647" "1648" "1649" "1650"
## [1651] "1651" "1652" "1653" "1654" "1655" "1656" "1657" "1658" "1659" "1660"
## [1661] "1661" "1662" "1663" "1664" "1665" "1666" "1667" "1668" "1669" "1670"
## [1671] "1671" "1672" "1673" "1674" "1675" "1676" "1677" "1678" "1679" "1680"
## [1681] "1681" "1682" "1683" "1684" "1685" "1686" "1687" "1688" "1689" "1690"
## [1691] "1691" "1692" "1693" "1694" "1695" "1696" "1697" "1698" "1699" "1700"
## [1701] "1701" "1702" "1703" "1704" "1705" "1706" "1707" "1708" "1709" "1710"
## [1711] "1711" "1712" "1713" "1714" "1715" "1716" "1717" "1718" "1719" "1720"
## [1721] "1721" "1722" "1723" "1724" "1725" "1726" "1727" "1728" "1729" "1730"
## [1731] "1731" "1732" "1733" "1734" "1735" "1736" "1737" "1738" "1739" "1740"
## [1741] "1741" "1742" "1743" "1744" "1745" "1746" "1747" "1748" "1749" "1750"
## [1751] "1751" "1752" "1753" "1754" "1755" "1756" "1757" "1758" "1759" "1760"
## [1761] "1761" "1762" "1763" "1764" "1765" "1766" "1767" "1768" "1769" "1770"
## [1771] "1771" "1772" "1773" "1774" "1775" "1776" "1777" "1778" "1779" "1780"
## [1781] "1781" "1782" "1783" "1784" "1785" "1786" "1787" "1788" "1789" "1790"
## [1791] "1791" "1792" "1793" "1794" "1795" "1796" "1797" "1798" "1799" "1800"
## [1801] "1801" "1802" "1803" "1804" "1805" "1806" "1807" "1808" "1809" "1810"
## [1811] "1811" "1812" "1813" "1814" "1815" "1816" "1817" "1818" "1819" "1820"
## [1821] "1821" "1822" "1823" "1824" "1825" "1826" "1827" "1828" "1829" "1830"
## [1831] "1831" "1832" "1833" "1834" "1835" "1836" "1837" "1838" "1839" "1840"
## [1841] "1841" "1842" "1843" "1844" "1845" "1846" "1847" "1848" "1849" "1850"
## [1851] "1851" "1852" "1853" "1854" "1855" "1856" "1857" "1858" "1859" "1860"
## [1861] "1861" "1862" "1863" "1864" "1865" "1866" "1867" "1868" "1869" "1870"
## [1871] "1871" "1872" "1873" "1874" "1875" "1876" "1877" "1878" "1879" "1880"
## [1881] "1881" "1882" "1883" "1884" "1885" "1886" "1887" "1888" "1889" "1890"
## [1891] "1891" "1892" "1893" "1894" "1895" "1896" "1897" "1898" "1899" "1900"
## [1901] "1901" "1902" "1903" "1904" "1905" "1906" "1907" "1908" "1909" "1910"
## [1911] "1911" "1912" "1913" "1914" "1915" "1916" "1917" "1918" "1919" "1920"
## [1921] "1921" "1922" "1923" "1924" "1925" "1926" "1927" "1928" "1929" "1930"
## [1931] "1931" "1932" "1933" "1934" "1935" "1936" "1937" "1938" "1939" "1940"
## [1941] "1941" "1942" "1943" "1944" "1945" "1946" "1947" "1948" "1949" "1950"
## [1951] "1951" "1952" "1953" "1954" "1955" "1956" "1957" "1958" "1959" "1960"
## [1961] "1961" "1962" "1963" "1964" "1965" "1966" "1967" "1968" "1969" "1970"
## [1971] "1971" "1972" "1973" "1974" "1975" "1976" "1977" "1978" "1979" "1980"
## [1981] "1981" "1982" "1983" "1984" "1985" "1986" "1987" "1988" "1989" "1990"
## [1991] "1991" "1992" "1993" "1994" "1995" "1996" "1997" "1998" "1999" "2000"
## [2001] "2001" "2002" "2003" "2004" "2005" "2006" "2007" "2008" "2009" "2010"
## [2011] "2011" "2012" "2013" "2014" "2015" "2016" "2017" "2018" "2019" "2020"
## [2021] "2021" "2022" "2023" "2024" "2025" "2026" "2027" "2028" "2029" "2030"
## [2031] "2031" "2032" "2033" "2034" "2035" "2036" "2037" "2038" "2039" "2040"
## [2041] "2041" "2042" "2043" "2044" "2045" "2046" "2047" "2048" "2049" "2050"
## [2051] "2051" "2052" "2053" "2054" "2055" "2056" "2057" "2058" "2059" "2060"
## [2061] "2061" "2062" "2063" "2064" "2065" "2066" "2067" "2068" "2069" "2070"
## [2071] "2071" "2072" "2073" "2074" "2075" "2076" "2077" "2078" "2079" "2080"
## [2081] "2081" "2082" "2083" "2084" "2085" "2086" "2087" "2088" "2089" "2090"
## [2091] "2091" "2092" "2093" "2094" "2095" "2096" "2097" "2098" "2099" "2100"
## [2101] "2101" "2102" "2103" "2104" "2105" "2106" "2107" "2108" "2109" "2110"
## [2111] "2111" "2112" "2113" "2114" "2115" "2116" "2117" "2118" "2119" "2120"
## [2121] "2121" "2122" "2123" "2124" "2125" "2126" "2127" "2128" "2129" "2130"
## [2131] "2131" "2132" "2133" "2134" "2135" "2136" "2137" "2138" "2139" "2140"
## [2141] "2141" "2142" "2143" "2144" "2145" "2146" "2147" "2148" "2149" "2150"
## [2151] "2151" "2152" "2153" "2154" "2155" "2156" "2157" "2158" "2159" "2160"
## [2161] "2161" "2162" "2163" "2164" "2165" "2166" "2167" "2168" "2169" "2170"
## [2171] "2171" "2172" "2173" "2174" "2175" "2176" "2177" "2178" "2179" "2180"
## [2181] "2181" "2182" "2183" "2184" "2185" "2186" "2187" "2188" "2189" "2190"
## [2191] "2191" "2192" "2193" "2194" "2195" "2196" "2197" "2198" "2199" "2200"
## [2201] "2201" "2202" "2203" "2204" "2205" "2206" "2207" "2208" "2209" "2210"
## [2211] "2211" "2212" "2213" "2214" "2215" "2216" "2217" "2218" "2219" "2220"
## [2221] "2221" "2222" "2223" "2224" "2225" "2226" "2227" "2228" "2229" "2230"
## [2231] "2231" "2232" "2233" "2234" "2235" "2236" "2237" "2238" "2239" "2240"
## [2241] "2241" "2242" "2243" "2244" "2245" "2246" "2247" "2248" "2249" "2250"
## [2251] "2251" "2252" "2253" "2254" "2255" "2256" "2257" "2258" "2259" "2260"
## [2261] "2261" "2262" "2263" "2264" "2265" "2266" "2267" "2268" "2269" "2270"
## [2271] "2271" "2272" "2273" "2274" "2275" "2276" "2277" "2278" "2279" "2280"
## [2281] "2281" "2282" "2283" "2284" "2285" "2286" "2287" "2288" "2289" "2290"
## [2291] "2291" "2292" "2293" "2294" "2295" "2296" "2297" "2298" "2299" "2300"
## [2301] "2301" "2302" "2303" "2304" "2305" "2306" "2307" "2308" "2309" "2310"
## [2311] "2311" "2312" "2313" "2314" "2315" "2316" "2317" "2318" "2319" "2320"
## [2321] "2321" "2322" "2323" "2324" "2325" "2326" "2327" "2328" "2329" "2330"
## [2331] "2331" "2332" "2333" "2334" "2335" "2336" "2337" "2338" "2339" "2340"
## [2341] "2341" "2342" "2343" "2344" "2345" "2346" "2347" "2348" "2349" "2350"
## [2351] "2351" "2352" "2353" "2354" "2355" "2356" "2357" "2358" "2359" "2360"
## [2361] "2361" "2362" "2363" "2364" "2365" "2366" "2367" "2368" "2369" "2370"
## [2371] "2371" "2372" "2373" "2374" "2375" "2376" "2377" "2378" "2379" "2380"
## [2381] "2381" "2382" "2383" "2384" "2385" "2386" "2387" "2388" "2389" "2390"
## [2391] "2391" "2392" "2393" "2394" "2395" "2396" "2397" "2398" "2399" "2400"
## [2401] "2401" "2402" "2403" "2404" "2405" "2406" "2407" "2408" "2409" "2410"
## [2411] "2411" "2412" "2413" "2414" "2415" "2416" "2417" "2418" "2419" "2420"
## [2421] "2421" "2422" "2423" "2424" "2425" "2426" "2427" "2428" "2429" "2430"
## [2431] "2431" "2432" "2433" "2434" "2435" "2436" "2437" "2438" "2439" "2440"
## [2441] "2441" "2442" "2443" "2444" "2445" "2446" "2447" "2448" "2449" "2450"
## [2451] "2451" "2452" "2453" "2454" "2455" "2456" "2457" "2458" "2459" "2460"
## [2461] "2461" "2462" "2463" "2464" "2465" "2466" "2467" "2468" "2469" "2470"
## [2471] "2471" "2472" "2473" "2474" "2475" "2476" "2477" "2478" "2479" "2480"
## [2481] "2481" "2482" "2483" "2484" "2485" "2486" "2487" "2488" "2489" "2490"
## [2491] "2491" "2492" "2493" "2494" "2495" "2496" "2497" "2498" "2499" "2500"
## [2501] "2501" "2502" "2503" "2504" "2505" "2506" "2507" "2508" "2509" "2510"
## [2511] "2511" "2512" "2513" "2514" "2515" "2516" "2517" "2518" "2519" "2520"
## [2521] "2521" "2522" "2523" "2524" "2525" "2526" "2527" "2528" "2529" "2530"
## [2531] "2531" "2532" "2533" "2534" "2535" "2536" "2537" "2538" "2539" "2540"
## [2541] "2541" "2542" "2543" "2544" "2545" "2546" "2547" "2548" "2549" "2550"
## [2551] "2551" "2552" "2553" "2554" "2555" "2556" "2557" "2558" "2559" "2560"
## [2561] "2561" "2562" "2563" "2564" "2565" "2566" "2567" "2568" "2569" "2570"
## [2571] "2571" "2572" "2573" "2574" "2575" "2576" "2577" "2578" "2579" "2580"
## [2581] "2581" "2582" "2583" "2584" "2585" "2586" "2587" "2588" "2589" "2590"
## [2591] "2591" "2592" "2593" "2594" "2595" "2596" "2597" "2598" "2599" "2600"
## [2601] "2601" "2602" "2603" "2604" "2605" "2606" "2607" "2608" "2609" "2610"
## [2611] "2611" "2612" "2613" "2614" "2615" "2616" "2617" "2618" "2619" "2620"
## [2621] "2621" "2622" "2623" "2624" "2625" "2626" "2627" "2628" "2629" "2630"
## [2631] "2631" "2632" "2633" "2634" "2635" "2636" "2637" "2638" "2639" "2640"
## [2641] "2641" "2642" "2643" "2644" "2645" "2646" "2647" "2648" "2649" "2650"
## [2651] "2651" "2652" "2653" "2654" "2655" "2656" "2657" "2658" "2659" "2660"
## [2661] "2661" "2662" "2663" "2664" "2665" "2666" "2667" "2668" "2669" "2670"
## [2671] "2671" "2672" "2673" "2674" "2675" "2676" "2677" "2678" "2679" "2680"
## [2681] "2681" "2682" "2683" "2684" "2685" "2686" "2687" "2688" "2689" "2690"
## [2691] "2691" "2692" "2693" "2694" "2695" "2696" "2697" "2698" "2699" "2700"
## [2701] "2701" "2702" "2703" "2704" "2705" "2706" "2707" "2708" "2709" "2710"
## [2711] "2711" "2712" "2713" "2714" "2715" "2716" "2717" "2718" "2719" "2720"
## [2721] "2721" "2722" "2723" "2724" "2725" "2726" "2727" "2728" "2729" "2730"
## [2731] "2731" "2732" "2733" "2734" "2735" "2736" "2737" "2738" "2739" "2740"
## [2741] "2741" "2742" "2743" "2744" "2745" "2746" "2747" "2748" "2749" "2750"
## [2751] "2751" "2752" "2753" "2754" "2755" "2756" "2757" "2758" "2759" "2760"
## [2761] "2761" "2762" "2763" "2764" "2765" "2766" "2767" "2768" "2769" "2770"
## [2771] "2771" "2772" "2773" "2774" "2775" "2776" "2777" "2778" "2779" "2780"
## [2781] "2781" "2782" "2783" "2784" "2785" "2786" "2787" "2788" "2789" "2790"
## [2791] "2791" "2792" "2793" "2794" "2795" "2796" "2797" "2798" "2799" "2800"
## [2801] "2801" "2802" "2803" "2804" "2805" "2806" "2807" "2808" "2809" "2810"
## [2811] "2811" "2812" "2813" "2814" "2815" "2816" "2817" "2818" "2819" "2820"
## [2821] "2821" "2822" "2823" "2824" "2825" "2826" "2827" "2828" "2829" "2830"
## [2831] "2831" "2832" "2833" "2834" "2835" "2836" "2837" "2838" "2839" "2840"
## [2841] "2841" "2842" "2843" "2844" "2845" "2846" "2847" "2848" "2849" "2850"
## [2851] "2851" "2852" "2853" "2854" "2855" "2856" "2857" "2858" "2859" "2860"
## [2861] "2861" "2862" "2863" "2864" "2865" "2866" "2867" "2868" "2869" "2870"
## [2871] "2871" "2872" "2873" "2874" "2875" "2876" "2877" "2878" "2879" "2880"
## [2881] "2881" "2882" "2883" "2884" "2885" "2886" "2887" "2888" "2889" "2890"
## [2891] "2891" "2892" "2893" "2894" "2895" "2896" "2897" "2898" "2899" "2900"
## [2901] "2901" "2902" "2903" "2904" "2905" "2906" "2907" "2908" "2909" "2910"
## [2911] "2911" "2912" "2913" "2914" "2915" "2916" "2917" "2918" "2919" "2920"
## [2921] "2921" "2922" "2923" "2924" "2925" "2926" "2927" "2928" "2929" "2930"
## [2931] "2931" "2932" "2933" "2934" "2935" "2936" "2937" "2938" "2939" "2940"
## [2941] "2941" "2942" "2943" "2944" "2945" "2946" "2947" "2948" "2949" "2950"
## [2951] "2951" "2952" "2953" "2954" "2955" "2956" "2957" "2958" "2959" "2960"
## [2961] "2961" "2962" "2963" "2964" "2965" "2966" "2967" "2968" "2969" "2970"
## [2971] "2971" "2972" "2973" "2974" "2975" "2976" "2977" "2978" "2979" "2980"
## [2981] "2981" "2982" "2983" "2984" "2985" "2986" "2987" "2988" "2989" "2990"
## [2991] "2991" "2992" "2993" "2994" "2995" "2996" "2997" "2998" "2999" "3000"
## [3001] "3001" "3002" "3003" "3004" "3005" "3006" "3007" "3008" "3009" "3010"
## [3011] "3011" "3012" "3013" "3014" "3015" "3016" "3017" "3018" "3019" "3020"
## [3021] "3021" "3022" "3023" "3024" "3025" "3026" "3027" "3028" "3029" "3030"
## [3031] "3031" "3032" "3033" "3034" "3035" "3036" "3037" "3038" "3039" "3040"
## [3041] "3041" "3042" "3043" "3044" "3045" "3046" "3047" "3048" "3049" "3050"
## [3051] "3051" "3052" "3053" "3054" "3055" "3056" "3057" "3058" "3059" "3060"
## [3061] "3061" "3062" "3063" "3064" "3065" "3066" "3067" "3068" "3069" "3070"
## [3071] "3071" "3072" "3073" "3074" "3075" "3076" "3077" "3078" "3079" "3080"
## [3081] "3081" "3082" "3083" "3084" "3085" "3086" "3087" "3088" "3089" "3090"
## [3091] "3091" "3092" "3093" "3094" "3095" "3096" "3097" "3098" "3099" "3100"
## [3101] "3101" "3102" "3103" "3104" "3105" "3106" "3107" "3108" "3109" "3110"
## [3111] "3111" "3112" "3113" "3114" "3115" "3116" "3117" "3118" "3119" "3120"
## [3121] "3121" "3122" "3123" "3124" "3125" "3126" "3127" "3128" "3129" "3130"
## [3131] "3131" "3132" "3133" "3134" "3135" "3136" "3137" "3138" "3139" "3140"
## [3141] "3141" "3142" "3143" "3144" "3145" "3146" "3147" "3148" "3149" "3150"
## [3151] "3151" "3152" "3153" "3154" "3155" "3156" "3157" "3158" "3159" "3160"
## [3161] "3161" "3162" "3163" "3164" "3165" "3166" "3167" "3168" "3169" "3170"
## [3171] "3171" "3172" "3173" "3174" "3175" "3176" "3177" "3178" "3179" "3180"
## [3181] "3181" "3182" "3183" "3184" "3185" "3186" "3187" "3188" "3189" "3190"
## [3191] "3191" "3192" "3193" "3194" "3195" "3196" "3197" "3198" "3199" "3200"
## [3201] "3201" "3202" "3203" "3204" "3205" "3206" "3207" "3208" "3209" "3210"
## [3211] "3211" "3212" "3213" "3214" "3215" "3216" "3217" "3218" "3219" "3220"
## [3221] "3221" "3222" "3223" "3224" "3225" "3226" "3227" "3228" "3229" "3230"
## [3231] "3231" "3232" "3233" "3234" "3235" "3236" "3237" "3238" "3239" "3240"
## [3241] "3241" "3242" "3243" "3244" "3245" "3246" "3247" "3248" "3249" "3250"
## [3251] "3251" "3252" "3253" "3254" "3255" "3256" "3257" "3258" "3259" "3260"
## [3261] "3261" "3262" "3263" "3264" "3265" "3266" "3267" "3268" "3269" "3270"
## [3271] "3271" "3272" "3273" "3274" "3275" "3276" "3277" "3278" "3279" "3280"
## [3281] "3281" "3282" "3283" "3284" "3285" "3286" "3287" "3288" "3289" "3290"
## [3291] "3291" "3292" "3293" "3294" "3295" "3296" "3297" "3298" "3299" "3300"
## [3301] "3301" "3302" "3303" "3304" "3305" "3306" "3307" "3308" "3309" "3310"
## [3311] "3311" "3312" "3313" "3314" "3315" "3316" "3317" "3318" "3319" "3320"
## [3321] "3321" "3322" "3323" "3324" "3325" "3326" "3327" "3328" "3329" "3330"
## [3331] "3331" "3332" "3333" "3334" "3335" "3336" "3337" "3338" "3339" "3340"
## [3341] "3341" "3342" "3343" "3344" "3345" "3346" "3347" "3348" "3349" "3350"
## [3351] "3351" "3352" "3353" "3354" "3355" "3356" "3357" "3358" "3359" "3360"
## [3361] "3361" "3362" "3363" "3364" "3365" "3366" "3367" "3368" "3369" "3370"
## [3371] "3371" "3372" "3373" "3374" "3375" "3376" "3377" "3378" "3379" "3380"
## [3381] "3381" "3382" "3383" "3384" "3385" "3386" "3387" "3388" "3389" "3390"
## [3391] "3391" "3392" "3393" "3394" "3395" "3396" "3397" "3398" "3399" "3400"
## [3401] "3401" "3402" "3403" "3404" "3405" "3406" "3407" "3408" "3409" "3410"
## [3411] "3411" "3412" "3413" "3414" "3415" "3416" "3417" "3418" "3419" "3420"
## [3421] "3421" "3422" "3423" "3424" "3425" "3426" "3427" "3428" "3429" "3430"
## [3431] "3431" "3432" "3433" "3434" "3435" "3436" "3437" "3438" "3439" "3440"
## [3441] "3441" "3442" "3443" "3444" "3445" "3446" "3447" "3448" "3449" "3450"
## [3451] "3451" "3452" "3453" "3454" "3455" "3456" "3457" "3458" "3459" "3460"
## [3461] "3461" "3462" "3463" "3464" "3465" "3466" "3467" "3468" "3469" "3470"
## [3471] "3471" "3472" "3473" "3474" "3475" "3476" "3477" "3478" "3479" "3480"
## [3481] "3481" "3482" "3483" "3484" "3485" "3486" "3487" "3488" "3489" "3490"
## [3491] "3491" "3492" "3493" "3494" "3495" "3496" "3497" "3498" "3499" "3500"
## [3501] "3501" "3502" "3503" "3504" "3505" "3506" "3507" "3508" "3509" "3510"
## [3511] "3511" "3512" "3513" "3514" "3515" "3516" "3517" "3518" "3519" "3520"
## [3521] "3521" "3522" "3523" "3524" "3525" "3526" "3527" "3528" "3529" "3530"
## [3531] "3531" "3532" "3533" "3534" "3535" "3536" "3537" "3538" "3539" "3540"
## [3541] "3541" "3542" "3543" "3544" "3545" "3546" "3547" "3548" "3549" "3550"
## [3551] "3551" "3552" "3553" "3554" "3555" "3556" "3557" "3558" "3559" "3560"
## [3561] "3561" "3562" "3563" "3564" "3565" "3566" "3567" "3568" "3569" "3570"
## [3571] "3571" "3572" "3573" "3574" "3575" "3576" "3577" "3578" "3579" "3580"
## [3581] "3581" "3582" "3583" "3584" "3585" "3586" "3587" "3588" "3589" "3590"
## [3591] "3591" "3592" "3593" "3594" "3595" "3596" "3597" "3598" "3599" "3600"
## [3601] "3601" "3602" "3603" "3604" "3605" "3606" "3607" "3608" "3609" "3610"
## [3611] "3611" "3612" "3613" "3614" "3615" "3616" "3617" "3618" "3619" "3620"
## [3621] "3621" "3622" "3623" "3624" "3625" "3626" "3627" "3628" "3629" "3630"
## [3631] "3631" "3632" "3633" "3634" "3635" "3636" "3637" "3638" "3639" "3640"
## [3641] "3641" "3642" "3643" "3644" "3645" "3646" "3647" "3648" "3649" "3650"
## [3651] "3651" "3652" "3653" "3654" "3655" "3656" "3657" "3658" "3659" "3660"
## [3661] "3661" "3662" "3663" "3664" "3665" "3666" "3667" "3668" "3669" "3670"
## [3671] "3671" "3672" "3673" "3674" "3675" "3676" "3677" "3678" "3679" "3680"
## [3681] "3681" "3682" "3683" "3684" "3685" "3686" "3687" "3688" "3689" "3690"
## [3691] "3691" "3692" "3693" "3694" "3695" "3696" "3697" "3698" "3699" "3700"
## [3701] "3701" "3702" "3703" "3704" "3705" "3706" "3707" "3708" "3709" "3710"
## [3711] "3711" "3712" "3713" "3714" "3715" "3716" "3717" "3718" "3719" "3720"
## [3721] "3721" "3722" "3723" "3724" "3725" "3726" "3727" "3728" "3729" "3730"
## [3731] "3731" "3732" "3733" "3734" "3735" "3736" "3737" "3738" "3739" "3740"
## [3741] "3741" "3742" "3743" "3744" "3745" "3746" "3747" "3748" "3749" "3750"
## [3751] "3751" "3752" "3753" "3754" "3755" "3756" "3757" "3758" "3759" "3760"
## [3761] "3761" "3762" "3763" "3764" "3765" "3766" "3767" "3768" "3769" "3770"
## [3771] "3771" "3772" "3773" "3774" "3775" "3776" "3777" "3778" "3779" "3780"
## [3781] "3781" "3782" "3783" "3784" "3785" "3786" "3787" "3788" "3789" "3790"
## [3791] "3791" "3792" "3793" "3794" "3795" "3796" "3797" "3798" "3799" "3800"
## [3801] "3801" "3802" "3803" "3804" "3805" "3806" "3807" "3808" "3809" "3810"
## [3811] "3811" "3812" "3813" "3814" "3815" "3816" "3817" "3818" "3819" "3820"
## [3821] "3821" "3822" "3823" "3824" "3825" "3826" "3827" "3828" "3829" "3830"
## [3831] "3831" "3832" "3833" "3834" "3835" "3836" "3837" "3838" "3839" "3840"
## [3841] "3841" "3842" "3843" "3844" "3845" "3846" "3847" "3848" "3849" "3850"
## [3851] "3851" "3852" "3853" "3854" "3855" "3856" "3857" "3858" "3859" "3860"
## [3861] "3861" "3862" "3863" "3864" "3865" "3866" "3867" "3868" "3869" "3870"
## [3871] "3871" "3872" "3873" "3874" "3875" "3876" "3877" "3878" "3879" "3880"
## [3881] "3881" "3882" "3883" "3884" "3885" "3886" "3887" "3888" "3889" "3890"
## [3891] "3891" "3892" "3893" "3894" "3895" "3896" "3897" "3898" "3899" "3900"
## [3901] "3901" "3902" "3903" "3904" "3905" "3906" "3907" "3908" "3909" "3910"
## [3911] "3911" "3912" "3913" "3914" "3915" "3916" "3917" "3918" "3919" "3920"
## [3921] "3921" "3922" "3923" "3924" "3925" "3926" "3927" "3928" "3929" "3930"
## [3931] "3931" "3932" "3933" "3934" "3935" "3936" "3937" "3938" "3939" "3940"
## [3941] "3941" "3942" "3943" "3944" "3945" "3946" "3947" "3948" "3949" "3950"
## [3951] "3951" "3952" "3953" "3954" "3955" "3956" "3957" "3958" "3959" "3960"
## [3961] "3961" "3962" "3963" "3964" "3965" "3966" "3967" "3968" "3969" "3970"
## [3971] "3971" "3972" "3973" "3974" "3975" "3976" "3977" "3978" "3979" "3980"
## [3981] "3981" "3982" "3983" "3984" "3985" "3986" "3987" "3988" "3989" "3990"
## [3991] "3991" "3992" "3993" "3994" "3995" "3996" "3997" "3998" "3999" "4000"
## [4001] "4001" "4002" "4003" "4004" "4005" "4006" "4007" "4008" "4009" "4010"
## [4011] "4011" "4012" "4013" "4014" "4015" "4016" "4017" "4018" "4019" "4020"
## [4021] "4021" "4022" "4023" "4024" "4025" "4026" "4027" "4028" "4029" "4030"
## [4031] "4031" "4032" "4033" "4034" "4035" "4036" "4037" "4038" "4039" "4040"
## [4041] "4041" "4042" "4043" "4044" "4045" "4046" "4047" "4048" "4049" "4050"
## [4051] "4051" "4052" "4053" "4054" "4055" "4056" "4057" "4058" "4059" "4060"
## [4061] "4061" "4062" "4063" "4064" "4065" "4066" "4067" "4068" "4069" "4070"
## [4071] "4071" "4072" "4073" "4074" "4075" "4076" "4077" "4078" "4079" "4080"
## [4081] "4081" "4082" "4083" "4084" "4085" "4086" "4087" "4088" "4089" "4090"
## [4091] "4091" "4092" "4093" "4094" "4095" "4096" "4097" "4098" "4099" "4100"
## [4101] "4101" "4102" "4103" "4104" "4105" "4106" "4107" "4108" "4109" "4110"
## [4111] "4111" "4112" "4113" "4114" "4115" "4116" "4117" "4118" "4119" "4120"
## [4121] "4121" "4122" "4123" "4124" "4125" "4126" "4127" "4128" "4129" "4130"
## [4131] "4131" "4132" "4133" "4134" "4135" "4136" "4137" "4138" "4139" "4140"
## [4141] "4141" "4142" "4143" "4144" "4145" "4146" "4147" "4148" "4149" "4150"
## [4151] "4151" "4152" "4153" "4154" "4155" "4156" "4157" "4158" "4159" "4160"
## [4161] "4161" "4162" "4163" "4164" "4165" "4166" "4167" "4168" "4169" "4170"
## [4171] "4171" "4172" "4173" "4174" "4175" "4176" "4177" "4178" "4179" "4180"
## [4181] "4181" "4182" "4183" "4184" "4185" "4186" "4187" "4188" "4189" "4190"
## [4191] "4191" "4192" "4193" "4194" "4195" "4196" "4197" "4198" "4199" "4200"
## [4201] "4201" "4202" "4203" "4204" "4205" "4206" "4207" "4208" "4209" "4210"
## [4211] "4211" "4212" "4213" "4214" "4215" "4216" "4217" "4218" "4219" "4220"
## [4221] "4221" "4222" "4223" "4224" "4225" "4226" "4227" "4228" "4229" "4230"
## [4231] "4231" "4232" "4233" "4234" "4235" "4236" "4237" "4238" "4239" "4240"
## [4241] "4241" "4242" "4243" "4244" "4245" "4246" "4247" "4248" "4249" "4250"
## [4251] "4251" "4252" "4253" "4254" "4255" "4256" "4257" "4258" "4259" "4260"
## [4261] "4261" "4262" "4263" "4264" "4265" "4266" "4267" "4268" "4269" "4270"
## [4271] "4271" "4272" "4273" "4274" "4275" "4276" "4277" "4278" "4279" "4280"
## [4281] "4281" "4282" "4283" "4284" "4285" "4286" "4287" "4288" "4289" "4290"
## [4291] "4291" "4292" "4293" "4294" "4295" "4296" "4297" "4298" "4299" "4300"
## [4301] "4301" "4302" "4303" "4304" "4305" "4306" "4307" "4308" "4309" "4310"
## [4311] "4311" "4312" "4313" "4314" "4315" "4316" "4317" "4318" "4319" "4320"
## [4321] "4321" "4322" "4323" "4324" "4325" "4326" "4327" "4328" "4329" "4330"
## [4331] "4331" "4332" "4333" "4334" "4335" "4336" "4337" "4338" "4339" "4340"
## [4341] "4341" "4342" "4343" "4344" "4345" "4346" "4347" "4348" "4349" "4350"
## [4351] "4351" "4352" "4353" "4354" "4355" "4356" "4357" "4358" "4359" "4360"
## [4361] "4361" "4362" "4363" "4364" "4365" "4366" "4367" "4368" "4369" "4370"
## [4371] "4371" "4372" "4373" "4374" "4375" "4376" "4377" "4378" "4379" "4380"
## [4381] "4381" "4382" "4383" "4384" "4385" "4386" "4387" "4388" "4389" "4390"
## [4391] "4391" "4392" "4393" "4394" "4395" "4396" "4397" "4398" "4399" "4400"
## [4401] "4401" "4402" "4403" "4404" "4405" "4406" "4407" "4408" "4409" "4410"
## [4411] "4411" "4412" "4413" "4414" "4415" "4416" "4417" "4418" "4419" "4420"
## [4421] "4421" "4422" "4423" "4424" "4425" "4426" "4427" "4428" "4429" "4430"
## [4431] "4431" "4432" "4433" "4434" "4435" "4436" "4437" "4438" "4439" "4440"
## [4441] "4441" "4442" "4443" "4444" "4445" "4446" "4447" "4448" "4449" "4450"
## [4451] "4451" "4452" "4453" "4454" "4455" "4456" "4457" "4458" "4459" "4460"
## [4461] "4461" "4462" "4463" "4464" "4465" "4466" "4467" "4468" "4469" "4470"
## [4471] "4471" "4472" "4473" "4474" "4475" "4476" "4477" "4478" "4479" "4480"
## [4481] "4481" "4482" "4483" "4484" "4485" "4486" "4487" "4488" "4489" "4490"
## [4491] "4491" "4492" "4493" "4494" "4495" "4496" "4497" "4498" "4499" "4500"
## [4501] "4501" "4502" "4503" "4504" "4505" "4506" "4507" "4508" "4509" "4510"
## [4511] "4511" "4512" "4513" "4514" "4515" "4516" "4517" "4518" "4519" "4520"
## [4521] "4521" "4522" "4523" "4524" "4525" "4526" "4527" "4528" "4529" "4530"
## [4531] "4531" "4532" "4533" "4534" "4535" "4536" "4537" "4538" "4539" "4540"
## [4541] "4541" "4542" "4543" "4544" "4545" "4546" "4547" "4548" "4549" "4550"
## [4551] "4551" "4552" "4553" "4554" "4555" "4556" "4557" "4558" "4559" "4560"
## [4561] "4561" "4562" "4563" "4564" "4565" "4566" "4567" "4568" "4569" "4570"
## [4571] "4571" "4572" "4573" "4574" "4575" "4576" "4577" "4578" "4579" "4580"
## [4581] "4581" "4582" "4583" "4584" "4585" "4586" "4587" "4588" "4589" "4590"
## [4591] "4591" "4592" "4593" "4594" "4595" "4596" "4597" "4598" "4599" "4600"
## [4601] "4601" "4602" "4603" "4604" "4605" "4606" "4607" "4608" "4609" "4610"
## [4611] "4611" "4612" "4613" "4614" "4615" "4616" "4617" "4618" "4619" "4620"
## [4621] "4621" "4622" "4623" "4624" "4625" "4626" "4627" "4628" "4629" "4630"
## [4631] "4631" "4632" "4633" "4634" "4635" "4636" "4637" "4638" "4639" "4640"
## [4641] "4641" "4642" "4643" "4644" "4645" "4646" "4647" "4648" "4649" "4650"
## [4651] "4651" "4652" "4653" "4654" "4655" "4656" "4657" "4658" "4659" "4660"
## [4661] "4661" "4662" "4663" "4664" "4665" "4666" "4667" "4668" "4669" "4670"
## [4671] "4671" "4672" "4673" "4674" "4675" "4676" "4677" "4678" "4679" "4680"
## [4681] "4681" "4682" "4683" "4684" "4685" "4686" "4687" "4688" "4689" "4690"
## [4691] "4691" "4692" "4693" "4694" "4695" "4696" "4697" "4698" "4699" "4700"
## [4701] "4701" "4702" "4703" "4704" "4705" "4706" "4707" "4708" "4709" "4710"
## [4711] "4711" "4712" "4713" "4714" "4715" "4716" "4717" "4718" "4719" "4720"
## [4721] "4721" "4722" "4723" "4724" "4725" "4726" "4727" "4728" "4729" "4730"
## [4731] "4731" "4732" "4733" "4734" "4735" "4736" "4737" "4738" "4739" "4740"
## [4741] "4741" "4742" "4743" "4744" "4745" "4746" "4747" "4748" "4749" "4750"
## [4751] "4751" "4752" "4753" "4754" "4755" "4756" "4757" "4758" "4759" "4760"
## [4761] "4761" "4762" "4763" "4764" "4765" "4766" "4767" "4768" "4769" "4770"
## [4771] "4771" "4772" "4773" "4774" "4775" "4776" "4777" "4778" "4779" "4780"
## [4781] "4781" "4782" "4783" "4784" "4785" "4786" "4787" "4788" "4789" "4790"
## [4791] "4791" "4792" "4793" "4794" "4795" "4796" "4797" "4798" "4799" "4800"
## [4801] "4801" "4802" "4803" "4804" "4805" "4806" "4807" "4808" "4809" "4810"
## [4811] "4811" "4812" "4813" "4814" "4815" "4816" "4817" "4818" "4819" "4820"
## [4821] "4821" "4822" "4823" "4824" "4825" "4826" "4827" "4828" "4829" "4830"
## [4831] "4831" "4832" "4833" "4834" "4835" "4836" "4837" "4838" "4839" "4840"
## [4841] "4841" "4842" "4843" "4844" "4845" "4846" "4847" "4848" "4849" "4850"
## [4851] "4851" "4852" "4853" "4854" "4855" "4856" "4857" "4858" "4859" "4860"
## [4861] "4861" "4862" "4863" "4864" "4865" "4866" "4867" "4868" "4869" "4870"
## [4871] "4871" "4872" "4873" "4874" "4875" "4876" "4877" "4878" "4879" "4880"
## [4881] "4881" "4882" "4883" "4884" "4885" "4886" "4887" "4888" "4889" "4890"
## [4891] "4891" "4892" "4893" "4894" "4895" "4896" "4897" "4898" "4899" "4900"
## [4901] "4901" "4902" "4903" "4904" "4905" "4906" "4907" "4908" "4909" "4910"
## [4911] "4911" "4912" "4913" "4914" "4915" "4916" "4917" "4918" "4919" "4920"
## [4921] "4921" "4922" "4923" "4924" "4925" "4926" "4927" "4928" "4929" "4930"
## [4931] "4931" "4932" "4933" "4934" "4935" "4936" "4937" "4938" "4939" "4940"
## [4941] "4941" "4942" "4943" "4944" "4945" "4946" "4947" "4948" "4949" "4950"
## [4951] "4951" "4952" "4953" "4954" "4955" "4956" "4957" "4958" "4959" "4960"
## [4961] "4961" "4962" "4963" "4964" "4965" "4966" "4967" "4968" "4969" "4970"
## [4971] "4971" "4972" "4973" "4974" "4975" "4976" "4977" "4978" "4979" "4980"
## [4981] "4981" "4982" "4983" "4984" "4985" "4986" "4987" "4988" "4989" "4990"
## [4991] "4991" "4992" "4993" "4994" "4995" "4996" "4997" "4998" "4999" "5000"
## [5001] "5001" "5002" "5003" "5004" "5005" "5006" "5007" "5008" "5009" "5010"
## [5011] "5011" "5012" "5013" "5014" "5015" "5016" "5017" "5018" "5019" "5020"
## [5021] "5021" "5022" "5023" "5024" "5025" "5026" "5027" "5028" "5029" "5030"
## [5031] "5031" "5032" "5033" "5034" "5035" "5036" "5037" "5038" "5039" "5040"
## [5041] "5041" "5042" "5043" "5044" "5045" "5046" "5047" "5048" "5049" "5050"
## [5051] "5051" "5052" "5053" "5054" "5055" "5056" "5057" "5058" "5059" "5060"
## [5061] "5061" "5062" "5063" "5064" "5065" "5066" "5067" "5068" "5069" "5070"
## [5071] "5071" "5072" "5073" "5074" "5075" "5076" "5077" "5078" "5079" "5080"
## [5081] "5081" "5082" "5083" "5084" "5085" "5086" "5087" "5088" "5089" "5090"
## [5091] "5091" "5092" "5093" "5094" "5095" "5096" "5097" "5098" "5099" "5100"
## [5101] "5101" "5102" "5103" "5104" "5105" "5106" "5107" "5108" "5109" "5110"
## [5111] "5111" "5112" "5113" "5114" "5115" "5116" "5117" "5118" "5119" "5120"
## [5121] "5121" "5122" "5123" "5124" "5125" "5126" "5127" "5128" "5129" "5130"
## [5131] "5131" "5132" "5133" "5134" "5135" "5136" "5137" "5138" "5139" "5140"
## [5141] "5141" "5142" "5143" "5144" "5145" "5146" "5147" "5148" "5149" "5150"
## [5151] "5151" "5152" "5153" "5154" "5155" "5156" "5157" "5158" "5159" "5160"
## [5161] "5161" "5162" "5163" "5164" "5165" "5166" "5167" "5168" "5169" "5170"
## [5171] "5171" "5172" "5173" "5174" "5175" "5176" "5177" "5178" "5179" "5180"
## [5181] "5181" "5182" "5183" "5184" "5185" "5186" "5187" "5188" "5189" "5190"
## [5191] "5191" "5192" "5193" "5194" "5195" "5196" "5197" "5198" "5199" "5200"
## [5201] "5201" "5202" "5203" "5204" "5205" "5206" "5207" "5208" "5209" "5210"
## [5211] "5211" "5212" "5213" "5214" "5215" "5216" "5217" "5218" "5219" "5220"
## [5221] "5221" "5222" "5223" "5224" "5225" "5226" "5227" "5228" "5229" "5230"
## [5231] "5231" "5232" "5233" "5234" "5235" "5236" "5237" "5238" "5239" "5240"
## [5241] "5241" "5242" "5243" "5244" "5245" "5246" "5247" "5248" "5249" "5250"
## [5251] "5251" "5252" "5253" "5254" "5255" "5256" "5257" "5258" "5259" "5260"
## [5261] "5261" "5262" "5263" "5264" "5265" "5266" "5267" "5268" "5269" "5270"
## [5271] "5271" "5272" "5273" "5274" "5275" "5276" "5277" "5278" "5279" "5280"
## [5281] "5281" "5282" "5283" "5284" "5285" "5286" "5287" "5288" "5289" "5290"
## [5291] "5291" "5292" "5293" "5294" "5295" "5296" "5297" "5298" "5299" "5300"
## [5301] "5301" "5302" "5303" "5304" "5305" "5306" "5307" "5308" "5309" "5310"
## [5311] "5311" "5312" "5313" "5314" "5315" "5316" "5317" "5318" "5319" "5320"
## [5321] "5321" "5322" "5323" "5324" "5325" "5326" "5327" "5328" "5329" "5330"
## [5331] "5331" "5332" "5333" "5334" "5335" "5336" "5337" "5338" "5339" "5340"
## [5341] "5341" "5342" "5343" "5344" "5345" "5346" "5347" "5348" "5349" "5350"
## [5351] "5351" "5352" "5353" "5354" "5355" "5356" "5357" "5358" "5359" "5360"
## [5361] "5361" "5362" "5363" "5364" "5365" "5366" "5367" "5368" "5369" "5370"
## [5371] "5371" "5372" "5373" "5374" "5375" "5376" "5377" "5378" "5379" "5380"
## [5381] "5381" "5382" "5383" "5384" "5385" "5386" "5387" "5388" "5389" "5390"
## [5391] "5391" "5392" "5393" "5394" "5395" "5396" "5397" "5398" "5399" "5400"
## [5401] "5401" "5402" "5403" "5404" "5405" "5406" "5407" "5408" "5409" "5410"
## [5411] "5411" "5412" "5413" "5414" "5415" "5416" "5417" "5418" "5419" "5420"
## [5421] "5421" "5422" "5423" "5424" "5425" "5426" "5427" "5428" "5429" "5430"
## [5431] "5431" "5432" "5433" "5434" "5435" "5436" "5437" "5438" "5439" "5440"
## [5441] "5441" "5442" "5443" "5444" "5445" "5446" "5447" "5448" "5449" "5450"
## [5451] "5451" "5452" "5453" "5454" "5455" "5456" "5457" "5458" "5459" "5460"
## [5461] "5461" "5462" "5463" "5464" "5465" "5466" "5467" "5468" "5469" "5470"
## [5471] "5471" "5472" "5473" "5474" "5475" "5476" "5477" "5478" "5479" "5480"
## [5481] "5481" "5482" "5483" "5484" "5485" "5486" "5487" "5488" "5489" "5490"
## [5491] "5491" "5492" "5493" "5494" "5495" "5496" "5497" "5498" "5499" "5500"
## [5501] "5501" "5502" "5503" "5504" "5505" "5506" "5507" "5508" "5509" "5510"
## [5511] "5511" "5512" "5513" "5514" "5515" "5516" "5517" "5518" "5519" "5520"
## [5521] "5521" "5522" "5523" "5524" "5525" "5526" "5527" "5528" "5529" "5530"
## [5531] "5531" "5532" "5533" "5534" "5535" "5536" "5537" "5538" "5539" "5540"
## [5541] "5541" "5542" "5543" "5544" "5545" "5546" "5547" "5548" "5549" "5550"
## [5551] "5551" "5552" "5553" "5554" "5555" "5556" "5557" "5558" "5559" "5560"
## [5561] "5561" "5562" "5563" "5564" "5565" "5566" "5567" "5568" "5569" "5570"
## [5571] "5571" "5572" "5573" "5574" "5575" "5576" "5577" "5578" "5579" "5580"
## [5581] "5581" "5582" "5583" "5584" "5585" "5586" "5587" "5588" "5589" "5590"
## [5591] "5591" "5592" "5593" "5594" "5595" "5596" "5597" "5598" "5599" "5600"
## [5601] "5601" "5602" "5603" "5604" "5605" "5606" "5607" "5608" "5609" "5610"
## [5611] "5611" "5612" "5613" "5614" "5615" "5616" "5617" "5618" "5619" "5620"
## [5621] "5621" "5622" "5623" "5624" "5625" "5626" "5627" "5628" "5629" "5630"
## [5631] "5631" "5632" "5633" "5634" "5635" "5636" "5637" "5638" "5639" "5640"
## [5641] "5641" "5642" "5643" "5644" "5645" "5646" "5647" "5648" "5649" "5650"
## [5651] "5651" "5652" "5653" "5654" "5655" "5656" "5657" "5658" "5659" "5660"
## [5661] "5661" "5662" "5663" "5664" "5665" "5666" "5667" "5668" "5669" "5670"
## [5671] "5671" "5672" "5673" "5674" "5675" "5676" "5677" "5678" "5679" "5680"
## [5681] "5681" "5682" "5683" "5684" "5685" "5686" "5687" "5688" "5689" "5690"
## [5691] "5691" "5692" "5693" "5694" "5695" "5696" "5697" "5698" "5699" "5700"
## [5701] "5701" "5702" "5703" "5704" "5705" "5706" "5707" "5708" "5709" "5710"
## [5711] "5711" "5712" "5713" "5714" "5715" "5716" "5717" "5718" "5719" "5720"
## [5721] "5721" "5722" "5723" "5724" "5725" "5726" "5727" "5728" "5729" "5730"
## [5731] "5731" "5732" "5733" "5734" "5735" "5736" "5737" "5738" "5739" "5740"
## [5741] "5741" "5742" "5743" "5744" "5745" "5746" "5747" "5748" "5749" "5750"
## [5751] "5751" "5752" "5753" "5754" "5755" "5756" "5757" "5758" "5759" "5760"
## [5761] "5761" "5762" "5763" "5764" "5765" "5766" "5767" "5768" "5769" "5770"
## [5771] "5771" "5772" "5773" "5774" "5775" "5776" "5777" "5778" "5779" "5780"
## [5781] "5781" "5782" "5783" "5784" "5785" "5786" "5787" "5788" "5789" "5790"
## [5791] "5791" "5792" "5793" "5794" "5795" "5796" "5797" "5798" "5799" "5800"
## [5801] "5801" "5802" "5803" "5804" "5805" "5806" "5807" "5808" "5809" "5810"
## [5811] "5811" "5812" "5813" "5814" "5815" "5816" "5817" "5818" "5819" "5820"
## [5821] "5821" "5822" "5823" "5824" "5825" "5826" "5827" "5828" "5829" "5830"
## [5831] "5831" "5832" "5833" "5834" "5835" "5836" "5837" "5838" "5839" "5840"
## [5841] "5841" "5842" "5843" "5844" "5845" "5846" "5847" "5848" "5849" "5850"
## [5851] "5851" "5852" "5853" "5854" "5855" "5856" "5857" "5858" "5859" "5860"
## [5861] "5861" "5862" "5863" "5864" "5865" "5866" "5867" "5868" "5869" "5870"
## [5871] "5871" "5872" "5873" "5874" "5875" "5876" "5877" "5878" "5879" "5880"
## [5881] "5881" "5882" "5883" "5884" "5885" "5886" "5887" "5888" "5889" "5890"
## [5891] "5891" "5892" "5893" "5894" "5895" "5896" "5897" "5898" "5899" "5900"
## [5901] "5901" "5902" "5903" "5904" "5905" "5906" "5907" "5908" "5909" "5910"
## [5911] "5911" "5912" "5913" "5914" "5915" "5916" "5917" "5918" "5919" "5920"
## [5921] "5921" "5922" "5923" "5924" "5925" "5926" "5927" "5928" "5929" "5930"
## [5931] "5931" "5932" "5933" "5934" "5935" "5936" "5937" "5938" "5939" "5940"
## [5941] "5941" "5942" "5943" "5944" "5945" "5946" "5947" "5948" "5949" "5950"
## [5951] "5951" "5952" "5953" "5954" "5955" "5956" "5957" "5958" "5959" "5960"
## [5961] "5961" "5962" "5963" "5964" "5965" "5966" "5967" "5968" "5969" "5970"
## [5971] "5971" "5972" "5973" "5974" "5975" "5976" "5977" "5978" "5979" "5980"
## [5981] "5981" "5982" "5983" "5984" "5985" "5986" "5987" "5988" "5989" "5990"
## [5991] "5991" "5992" "5993" "5994" "5995" "5996" "5997" "5998" "5999" "6000"
##
## $chain
## [1] "1" "2" "3" "4"
##
## $variable
## [1] "b_Intercept" "b_sigma_Intercept" "b_sFat"
## [4] "b_sAge" "b_breedckcs" "b_breedpug"
## [7] "b_breedret" "b_breedother" "b_sexMale"
## [10] "b_comorbidityyes" "b_sigma_comorbidityyes" "sd_visit__Intercept"
## [13] "alpha" "Intercept" "Intercept_sigma"
## [16] "r_visit[v0,Intercept]" "r_visit[v1,Intercept]" "prior_Intercept"
## [19] "prior_b_sFat" "prior_b_sAge" "prior_b_breedckcs"
## [22] "prior_b_breedpug" "prior_b_breedret" "prior_b_breedother"
## [25] "prior_b_sexMale" "prior_b_comorbidityyes" "prior_Intercept_sigma"
## [28] "prior_alpha" "prior_sd_visit" "lprior"
## [31] "lp__" "z_1[1,1]" "z_1[1,2]"
mcmc_violin(posterior, pars = "b_sFat")
color_scheme_set("viridis")
mcmc_trace(posterior, pars = "b_sFat") +
theme_classic() +
labs(title = "Trace plot for beta of standardised fat",
y = "b_sFat",
x = "Rank") +
theme(axis.title.y = element_text(size=12, face="bold"),
axis.title.x = element_text(size=12, face="bold"),
title = element_text(size=12, face="bold"),
plot.title = element_text(hjust = 0.5),
axis.text.x = element_text(color = "grey50", size = 12),
axis.text.y = element_text(color = "grey8",size = 12))
mcmc_plot(model, type = 'rank_overlay', variable = "b_sFat") +
theme_classic() +
labs(title = "Trace rank plot for beta of standardised fat",
y = "b_sFat",
x = "Rank") +
ylim(250, 350) +
theme(axis.title.y = element_text(size=12, face="bold"),
axis.title.x = element_text(size=12, face="bold"),
title = element_text(size=12, face="bold"),
plot.title = element_text(hjust = 0.5),
axis.text.x = element_text(color = "grey50", size = 12),
axis.text.y = element_text(color = "grey8",size = 12))
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
Usually better than the compatoability intervals given in the summary
draws <- as.matrix(model)
HPDI(draws[,3], 0.97) # 1st column is draws for bf
## |0.97 0.97|
## 0.01884829 0.58929953
bayes_R2(model, probs = c(0.015, 0.5, 0.985)) # 0.015, 0.5, 0.985 are the quantiles
## Estimate Est.Error Q1.5 Q50 Q98.5
## R2 0.1835293 0.04804434 0.07906287 0.1837382 0.2893523
loo_R2(model, probs = c(0.015, 0.5, 0.985)) # 0.015, 0.5, 0.985 are the quantiles
## Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.
## Estimate Est.Error Q1.5 Q50 Q98.5
## R2 -0.1560929 0.09813997 -0.3939027 -0.1521871 0.04130273
checks whether actual data is similar to simulated data. ### a. standard
color_scheme_set("blue")
pp_check(model, ndraws = 100) +
theme_classic() +
labs(title = "Posterior predictions of body fat model (skew-normal)",
y = "Density",
x = "Log hair cortisol (standardised") +
theme(axis.title.y = element_text(size=12, face="bold"),
axis.title.x = element_text(size=12, face="bold"),
title = element_text(size=12, face="bold"),
plot.title = element_text(hjust = 0.5),
axis.text.x = element_text(color = "grey50", size = 12),
axis.text.y = element_text(color = "grey8",size = 12), ,
legend.position = "inside", legend.position.inside = c(0.91, 0.83))
color_scheme_set("blue")
pp_check(model, ndraws = 100) +
theme_classic() +
labs(title = "Posterior predictions of body fat model (skew-normal)",
y = "Density",
x = "Log hair cortisol (standardised") +
xlim(-5, 7) +
theme(axis.title.y = element_text(size=12, face="bold"),
axis.title.x = element_text(size=12, face="bold"),
title = element_text(size=12, face="bold"),
plot.title = element_text(hjust = 0.5),
axis.text.x = element_text(color = "grey50", size = 12),
axis.text.y = element_text(color = "grey8",size = 12), ,
legend.position = "inside", legend.position.inside = c(0.91, 0.83))
par(mfrow = c(1,1))
pp_check(model, type = "hist", ndraws = 11, binwidth = 0.25) # separate histograms of 11 MCMC draws vs actual data
pp_check(model, type = "error_hist", ndraws = 11) # separate histograms of errors for 11 draws
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
pp_check(model, type = "scatter_avg", ndraws = 100) # scatter plot
pp_check(model, type = "stat_2d") # scatterplot of joint posteriors
## Using all posterior draws for ppc type 'stat_2d' by default.
## Note: in most cases the default test statistic 'mean' is too weak to detect anything of interest.
# PPC functions for predictive checks based on (approximate) leave-one-out (LOO) cross-validation
pp_check(model, type = "loo_pit_overlay", ndraws = 1000)
## Warning: Found 5 observations with a pareto_k > 0.67 in model '.x1'. We
## recommend to run more iterations to get at least about 2200 posterior draws to
## improve LOO-CV approximation accuracy.
## NOTE: The kernel density estimate assumes continuous observations and is not optimal for discrete observations.
# separate histograms of 11 MCMC draws vs actual data
pp_check(model, type = "hist", ndraws = 11, binwidth = 0.25) +
theme_classic() +
labs(title = "Actual data vs 11 MCMC draws for body fat model") +
theme(axis.title.y = element_text(size=12, face="bold"),
axis.title.x = element_text(size=12, face="bold"),
title = element_text(size=12, face="bold"),
plot.title = element_text(hjust = 0.5),
axis.text.x = element_text(color = "grey50", size = 12),
axis.text.y = element_text(color = "grey8",size = 12))
# scatterplot of joint posteriors
pp_check(model, type = "stat_2d") +
theme_classic() +
labs(title = "Scatterplot of joint posteriors for body fat model") +
theme(axis.title.y = element_text(size=12, face="bold"),
axis.title.x = element_text(size=12, face="bold"),
title = element_text(size=12, face="bold"),
plot.title = element_text(hjust = 0.5),
axis.text.x = element_text(color = "grey50", size = 12),
axis.text.y = element_text(color = "grey8",size = 12), ,
legend.position = "inside", legend.position.inside = c(0.91, 0.13))
## Using all posterior draws for ppc type 'stat_2d' by default.
## Note: in most cases the default test statistic 'mean' is too weak to detect anything of interest.
pp_check(model, type = "loo_pit_overlay", ndraws = 1000) +
theme_classic() +
labs(title = "'LOO PIT overlay' plot for body fat model") +
theme(axis.title.y = element_text(size=12, face="bold"),
axis.title.x = element_text(size=12, face="bold"),
title = element_text(size=12, face="bold"),
plot.title = element_text(hjust = 0.5),
axis.text.x = element_text(color = "grey50", size = 12),
axis.text.y = element_text(color = "grey8",size = 12),
legend.position = "inside", legend.position.inside = c(0.91, 0.13))
## Warning: Found 7 observations with a pareto_k > 0.67 in model '.x1'. We
## recommend to run more iterations to get at least about 2200 posterior draws to
## improve LOO-CV approximation accuracy.
## NOTE: The kernel density estimate assumes continuous observations and is not optimal for discrete observations.
pp_check(model, type = "error_scatter_avg")
## Using all posterior draws for ppc type 'error_scatter_avg' by default.
color_scheme_set("blue")
pairs(model, regex = TRUE)
color_scheme_set("pink")
mcmc_pairs(posterior,
pars = c("Intercept", "Intercept_sigma",
"alpha",
"b_sFat", "b_sAge",
"b_sexMale",
"b_comorbidityyes", "b_sigma_comorbidityyes"),
off_diag_args = list(size = 0.5))
color_scheme_set("pink")
mcmc_pairs(posterior,
pars = c("Intercept", "Intercept_sigma",
"alpha",
"b_breedckcs", "b_breedpug",
"b_breedret", "b_breedother"),
off_diag_args = list(size = 0.5))
loo_model <- loo(model, moment_match = TRUE)
loo_model
##
## Computed from 24000 by 55 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -80.4 5.5
## p_loo 11.3 2.4
## looic 160.8 11.0
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 1.2]).
##
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
First, check the sensitivity of the prior and likelihood to power-scaling. Posterior and posteriors resulting from power-scaling.
color_scheme_set("blue")
powerscale_sensitivity(model, variable = c("b_Intercept",
"b_sFat",
"b_sAge",
"b_breedckcs", "b_breedother", "b_breedpug", "b_breedret",
"b_sexMale",
"b_comorbidityyes"))
## Sensitivity based on cjs_dist
## Prior selection: all priors
## Likelihood selection: all data
##
## variable prior likelihood diagnosis
## b_Intercept 0.042 0.053 -
## b_sFat 0.032 0.213 -
## b_sAge 0.018 0.133 -
## b_breedckcs 0.027 0.070 -
## b_breedother 0.025 0.089 -
## b_breedpug 0.025 0.112 -
## b_breedret 0.027 0.092 -
## b_sexMale 0.015 0.094 -
## b_comorbidityyes 0.026 0.165 -
check_prior(model, effects = "all")
## Parameter Prior_Quality
## 1 b_Intercept informative
## 2 b_sFat informative
## 3 b_sAge informative
## 4 b_breedckcs informative
## 5 b_breedpug informative
## 6 b_breedret informative
## 7 b_breedother informative
## 8 b_sexMale informative
## 9 b_comorbidityyes informative
## 10 sd_visit__Intercept informative
prior <- prior_draws(model)
prior %>% glimpse()
## Rows: 24,000
## Columns: 12
## $ Intercept <dbl> 1.0290702041, 0.3306087896, 0.4839608205, 0.008573827…
## $ b_sFat <dbl> 0.03607712, -0.54894592, 0.77995266, 0.83958415, -0.0…
## $ b_sAge <dbl> -0.13587838, 0.21510924, 0.02527679, -0.81112013, 0.4…
## $ b_breedckcs <dbl> 0.47432534, -2.09018501, 0.81528623, 0.42376328, -0.6…
## $ b_breedpug <dbl> 1.79481210, -1.08212312, -2.04839099, 0.86380913, 0.5…
## $ b_breedret <dbl> 0.9778925, -1.1427470, 0.9106284, -0.3588809, -1.0509…
## $ b_breedother <dbl> -0.10398011, 1.11981942, 0.20578397, -2.64647193, 1.2…
## $ b_sexMale <dbl> -1.260058872, 1.730863388, 1.910090289, -2.017788670,…
## $ b_comorbidityyes <dbl> -0.332723946, -0.002110186, 0.912710899, -0.232086474…
## $ Intercept_sigma <dbl> -2.2490940, 3.9302854, -2.3171893, -5.9297446, 10.597…
## $ alpha <dbl> 4.43100344, 3.73963395, 5.94892021, 5.50938705, 4.706…
## $ sd_visit <dbl> 0.7810777, 0.1843761, 1.0139941, 1.1335228, 1.0756120…
set.seed(5)
prior %>%
slice_sample(n = 50) %>%
rownames_to_column("draw") %>%
expand_grid(a = c(-2, 2)) %>%
mutate(c = Intercept + b_sAge * a) %>%
ggplot(aes(x = a, y = c)) +
geom_line(aes(group = draw),
color = "firebrick", alpha = .4) +
labs(x = "Age (std)",
y = "log(cort) (std)") +
coord_cartesian(ylim = c(-4, 4)) +
theme_bw() +
theme(panel.grid = element_blank())
set.seed(5)
prior %>%
slice_sample(n = 50) %>%
rownames_to_column("draw") %>%
expand_grid(a = c(-2, 2)) %>%
mutate(c = Intercept + b_sFat * a) %>%
ggplot(aes(x = a, y = c)) +
geom_line(aes(group = draw),
color = "firebrick", alpha = .4) +
labs(x = "Fat (std)",
y = "log(cort) (std)") +
coord_cartesian(ylim = c(-2, 2)) +
theme_bw() +
theme(panel.grid = element_blank())
set.seed(5)
prior %>%
slice_sample(n = 50) %>%
rownames_to_column("draw") %>%
expand_grid(a = c(-2, 2)) %>%
mutate(c = Intercept + b_sFat * a) %>%
ggplot(aes(x = a, y = c)) +
geom_line(aes(group = draw),
color = "firebrick", alpha = .4) +
labs(x = "Fat (standardised)",
y = "Log hair cortisol (standardisedd)") +
coord_cartesian(ylim = c(-2, 2)) +
theme_classic() +
theme(panel.grid = element_blank()) +
labs(title = "Draws from prior for effect of body fat on hair cortisol") +
theme(axis.title.y = element_text(size=12, face="bold"),
axis.title.x = element_text(size=12, face="bold"),
title = element_text(size=12, face="bold"),
plot.title = element_text(hjust = 0.5),
axis.text.x = element_text(color = "grey50", size = 12),
axis.text.y = element_text(color = "grey8",size = 12))
set.seed(5)
post_draws <- draws
post_draws <- as.data.frame(post_draws)
post_draws %>%
slice_sample(n = 50) %>%
rownames_to_column("draw") %>%
expand_grid(a = c(-2, 2)) %>%
mutate(c = Intercept + b_sFat * a) %>%
ggplot(aes(x = a, y = c)) +
geom_line(aes(group = draw),
color = "steelblue3", alpha = .4) +
labs(x = "Fat (standardised)",
y = "Log hair cortisol (standardisedd)") +
coord_cartesian(ylim = c(-2, 2)) +
theme_classic() +
theme(panel.grid = element_blank()) +
labs(title = "Draws from posterior for effect of body fat on hair cortisol") +
theme(axis.title.y = element_text(size=12, face="bold"),
axis.title.x = element_text(size=12, face="bold"),
title = element_text(size=12, face="bold"),
plot.title = element_text(hjust = 0.5),
axis.text.x = element_text(color = "grey50", size = 12),
axis.text.y = element_text(color = "grey8",size = 12))
set.seed(5)
prior %>%
slice_sample(n = 50) %>%
rownames_to_column("draw") %>%
expand_grid(a = c(0, 1)) %>%
mutate(c = Intercept + b_comorbidityyes * a) %>%
ggplot(aes(x = a, y = c)) +
geom_line(aes(group = draw),
color = "firebrick", alpha = .4) +
geom_point(color = "firebrick", size = 2) +
labs(x = "comorbidity",
y = "log(cort) (std)") +
coord_cartesian(ylim = c(-2, 2)) +
theme_bw() +
theme(panel.grid = element_blank())
set.seed(5)
prior %>%
slice_sample(n = 50) %>%
rownames_to_column("draw") %>%
expand_grid(a = c(0, 1)) %>%
mutate(c = Intercept + b_sexMale * a) %>%
ggplot(aes(x = a, y = c)) +
geom_line(aes(group = draw),
color = "firebrick", alpha = .4) +
geom_point(color = "firebrick", size = 2) +
labs(x = "Sex breed",
y = "log(cort) (std)") +
coord_cartesian(ylim = c(-2, 2)) +
theme_bw() +
theme(panel.grid = element_blank())
set.seed(5)
prior %>%
slice_sample(n = 50) %>%
rownames_to_column("draw") %>%
expand_grid(a = c(0, 1)) %>%
mutate(c = Intercept + b_breedckcs * a) %>%
ggplot(aes(x = a, y = c)) +
geom_line(aes(group = draw),
color = "firebrick", alpha = .4) +
geom_point(color = "firebrick", size = 2) +
labs(x = "CKCS breed",
y = "log(cort) (std)") +
coord_cartesian(ylim = c(-2, 2)) +
theme_bw() +
theme(panel.grid = element_blank())
set.seed(5)
prior %>%
slice_sample(n = 50) %>%
rownames_to_column("draw") %>%
expand_grid(a = c(0, 1)) %>%
mutate(c = Intercept + b_breedother * a) %>%
ggplot(aes(x = a, y = c)) +
geom_line(aes(group = draw),
color = "firebrick", alpha = .4) +
geom_point(color = "firebrick", size = 2) +
labs(x = "Other breed",
y = "log(cort) (std)") +
coord_cartesian(ylim = c(-4, 4)) +
theme_bw() +
theme(panel.grid = element_blank())
set.seed(5)
prior %>%
slice_sample(n = 50) %>%
rownames_to_column("draw") %>%
expand_grid(a = c(0, 1)) %>%
mutate(c = Intercept + b_breedpug * a) %>%
ggplot(aes(x = a, y = c)) +
geom_line(aes(group = draw),
color = "firebrick", alpha = .4) +
geom_point(color = "firebrick", size = 2) +
labs(x = "Pug breed",
y = "log(cort) (std)") +
coord_cartesian(ylim = c(-4, 4)) +
theme_bw() +
theme(panel.grid = element_blank())
set.seed(5)
prior %>%
slice_sample(n = 50) %>%
rownames_to_column("draw") %>%
expand_grid(a = c(0, 1)) %>%
mutate(c = Intercept + b_breedpug * a) %>%
ggplot(aes(x = a, y = c)) +
geom_line(aes(group = draw),
color = "firebrick", alpha = .4) +
geom_point(color = "firebrick", size = 2) +
labs(x = "Pug breed",
y = "log(cort) (std)") +
coord_cartesian(ylim = c(-4, 4)) +
theme_bw() +
theme(panel.grid = element_blank())
Can simulate data just on the priors. Fit model but only consider prior when fitting model. If this looks reasonable, it helps to confirm that your priors were reasonable
set.seed(666)
model_priors_only <- brm(slgCort ~ sFat + sAge + breed + sex + comorbidity + (1 | visit),
family = skew_normal(),
prior = priors,
data = df2,
sample_prior = "only",
save_pars = save_pars(all =TRUE))
## Compiling Stan program...
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 17.0.0 (clang-1700.0.13.5)’
## using SDK: ‘MacOSX15.5.sdk’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/opt/R/arm64/include -fPIC -falign-functions=64 -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
## 679 | #include <cmath>
## | ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1
## Start sampling
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 7.6e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.76 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.034 seconds (Warm-up)
## Chain 1: 0.033 seconds (Sampling)
## Chain 1: 0.067 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 4e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.035 seconds (Warm-up)
## Chain 2: 0.028 seconds (Sampling)
## Chain 2: 0.063 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 5e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.05 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.038 seconds (Warm-up)
## Chain 3: 0.033 seconds (Sampling)
## Chain 3: 0.071 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 3e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.037 seconds (Warm-up)
## Chain 4: 0.027 seconds (Sampling)
## Chain 4: 0.064 seconds (Total)
## Chain 4:
## Warning: There were 1 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: Examine the pairs() plot to diagnose sampling problems
color_scheme_set("pink")
pp_check(model_priors_only, ndraws = 100) +
xlim(-10, 20) +
theme_classic() +
labs(title = "Prior predictions of body fat model versus data",
y = "Density",
x = "Log hair cortisol (standardised") +
theme(axis.title.y = element_text(size=12, face="bold"),
axis.title.x = element_text(size=12, face="bold"),
title = element_text(size=12, face="bold"),
plot.title = element_text(hjust = 0.5),
axis.text.x = element_text(color = "grey50", size = 12),
axis.text.y = element_text(color = "grey8",size = 12),
legend.position = "inside", legend.position.inside = c(0.91, 0.83))
## Warning: Removed 74 rows containing non-finite outside the scale range
## (`stat_density()`).
NB Uses posterior_predict ## 1. Posterior predictive distribution plots for a continuous predictor variable
par(mfrow = c(2,2))
# plot the observed data
plot(slgCort ~ sFat, main = "Observed", xlab = "Fat (std)", ylab = "log hair cortisol (std)",
data = df2, ylim = c(-3, 3), xlim = c(-3,3), col = "steelblue3") # This is the observed data
# use posterior predict to simulate predictions
ppd <- posterior_predict(model)
# Plot 3 simulations of the data
plot(ppd[50,] ~ df2$sFat, main = "Simulated", xlab = "Fat (std)", ylab = "log hair cortisol (std)",
ylim = c(-3,3), xlim = c(-3,3), col = "firebrick")
plot(ppd[51,] ~ df2$sFat, main = "Simulated", xlab = "Fat (std)", ylab = "log hair cortisol (std)",
ylim = c(-3,3), xlim = c(-3,3), col = "firebrick")
plot(ppd[52,] ~ df2$sFat, main = "Simulated", xlab = "Fat (std)", ylab = "log hair cortisol (std)",
ylim = c(-3,3), xlim = c(-3,3), col = "firebrick")
par(mfrow = c(2,2))
vioplot(slgCort ~ comorbidity, data = df2, main = "Observed", col = "steelblue")
vioplot(ppd[sample(seq(1, dim(ppd)[1]), 1),] ~ df2$comorbidity, main = "PPD", col = "firebrick")
vioplot(ppd[sample(seq(1, dim(ppd)[1]), 1),] ~ df2$comorbidity, main = "PPD", col = "firebrick")
vioplot(ppd[sample(seq(1, dim(ppd)[1]), 1),] ~ df2$comorbidity, main = "PPD", col = "firebrick")
par(mfrow = c(2,2))
vioplot(slgCort ~ sex, data = df2, main = "Observed", col = "steelblue3")
vioplot(ppd[sample(seq(1, dim(ppd)[1]), 1),] ~ df2$sex, main = "PPD", col = "firebrick")
vioplot(ppd[sample(seq(1, dim(ppd)[1]), 1),] ~ df2$sex, main = "PPD", col = "firebrick")
vioplot(ppd[sample(seq(1, dim(ppd)[1]), 1),] ~ df2$sex, main = "PPD", col = "firebrick")
par(mfrow = c(2,2))
stripchart(slgCort ~ breed, vertical = TRUE, method = "jitter",
col = "steelblue3", data = df2, pch = 20, main = "Observed")
stripchart(ppd[sample(seq(1, dim(ppd)[1]), 1),] ~ breed, vertical = TRUE, method = "jitter", col = "firebrick", data = df2, pch = 20, main = "PPD")
stripchart(ppd[sample(seq(1, dim(ppd)[1]), 1),] ~ breed, vertical = TRUE, method = "jitter", col = "firebrick", data = df2, pch = 20, main = "PPD")
stripchart(ppd[sample(seq(1, dim(ppd)[1]), 1),] ~ breed, vertical = TRUE, method = "jitter", col = "firebrick", data = df2, pch = 20, main = "PPD")
plot(conditional_effects(model), ask = FALSE)
ce <- conditional_effects(model, "sFat", method = "posterior_predict", prob = 0.95)
plot(ce, plot = FALSE)[[1]] +
theme_bw() +
labs(title = "Conditional effect of body fat on hair cortisol") +
labs(y = paste0("Log hair cortisol (standardised)")) +
labs(x = paste0("Body fat (standardised)")) +
theme(axis.title.y = element_text(size=12, face="bold"),
axis.title.x = element_text(size=12, face="bold"),
title = element_text(size=12, face="bold"),
plot.title = element_text(hjust = 0.5),
axis.text.x = element_text(color = "grey50", size = 12),
axis.text.y = element_text(color = "grey50", size = 12),
legend.position = "inside", legend.position.inside = c(0.93, 0.80))
ce <- conditional_effects(model, "sFat:comorbidity")
plot(ce, plot = FALSE)[[1]] +
theme_bw() +
labs(title = "Conditional effect of body fat on hair cortisol") +
labs(y = paste0("Log hair cortisol (standardised)")) +
labs(x = paste0("Body fat (standardised)")) +
theme(axis.title.y = element_text(size=12, face="bold"),
axis.title.x = element_text(size=12, face="bold"),
title = element_text(size=12, face="bold"),
plot.title = element_text(hjust = 0.5),
axis.text.x = element_text(color = "grey50", size = 12),
axis.text.y = element_text(color = "grey50", size = 12),
legend.position = "inside", legend.position.inside = c(0.91, 0.13))
color_scheme_set("blue")
mcmc_plot(model)
mcmc_plot(model,
variable = c("b_sFat", "prior_b_sFat"))
mcmc_plot(model,
variable = c("b_sFat", "prior_b_sFat"),
type = "areas") +
theme_classic() +
labs(title = "Prior vs posterior distribution for body fat effect") +
labs(y = "") +
labs(x = paste0("Possible parameter values")) +
scale_y_discrete(labels=c("prior_b_sFat" = "Prior", "b_sFat" = "Posterior"),
limits = c("prior_b_sFat", "b_sFat")) +
theme(axis.title.y = element_text(size=12, face="bold"),
axis.title.x = element_text(size=12, face="bold"),
title = element_text(size=12, face="bold"),
plot.title = element_text(hjust = 0.5),
axis.text.x = element_text(color = "grey50", size = 12),
axis.text.y = element_text(color = "grey8",size = 12))
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
mcmc_plot(model,
variable = "b_sexMale")
mcmc_plot(model,
variable = c("b_Intercept", "alpha",
"b_sFat",
"b_sAge",
"b_breedckcs",
"b_breedpug",
"b_breedret",
"b_breedother",
"b_sexMale",
"b_comorbidityyes"))
posterior <- as.matrix(model)
mcmc_areas(posterior,
pars = c("b_Intercept", "alpha",
"b_sFat",
"b_sAge",
"b_breedckcs",
"b_breedpug",
"b_breedret",
"b_breedother",
"b_sexMale",
"b_comorbidityyes"),
# arbitrary threshold for shading probability mass
prob = 0.75)
posterior <- as.matrix(model)
mcmc_areas(posterior,
pars = c("b_sFat",
"b_sAge",
"b_breedckcs",
"b_breedpug",
"b_breedret",
"b_breedother",
"b_sexMale",
"b_comorbidityyes"),
prob = 0.75) # arbitrary threshold for shading probability mass
posterior <- as.matrix(model)
mcmc_areas(posterior,
pars = "b_sFat",
# arbitrary threshold for shading probability mass
prob = 0.97) +
theme_classic() +
labs(title = "Posterior distribution for body fat effect",
y = "Density distribution",
x = "Possible parameter values") +
scale_y_discrete(labels=("b_sFat" = "Beta for fat mass")) +
theme(axis.title.y = element_text(size=12, face="bold"),
axis.title.x = element_text(size=12, face="bold"),
title = element_text(size=12, face="bold"),
plot.title = element_text(hjust = 0.5),
axis.text.x = element_text(color = "grey50", size = 12),
axis.text.y = element_text(color = "grey8",size = 12))
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
# Focus on describing posterior
hdi_range <- hdi(model, ci = c(0.65, 0.70, 0.80, 0.89, 0.95))
plot(hdi_range, show_intercept = T)
# Focus on describing posterior
hdi_range <- hdi(model, ci = c(0.65, 0.70, 0.80, 0.89, 0.97), parameters = "sFat")
plot(hdi_range, show_intercept = T) +
labs(title = "Posterior distribution for body fat effect") +
labs(y = "Density distribution") +
labs(x = "Possible parameter values") +
theme(axis.title.y = element_text(size=12, face="bold"),
axis.title.x = element_text(size=12, face="bold"),
title = element_text(size=12, face="bold"),
plot.title = element_text(hjust = 0.5),
axis.text.x = element_text(color = "grey50", size = 12),
axis.text.y = element_text(color = "grey8",size = 12))
##. a. Using ‘fitted’ which is a wrapper for posterior_epred (only returns error from model)
# Create new data for an average dog without a comorbidity
nd <- tibble(sFat = seq(from = -3.2, to = 3.2, length.out = 30), visit = "v0", breed = "mix", sex = "Female", comorbidity = "no", sAge = 0)
# now use `fitted()` to get the model-implied trajectories
pr1 <- fitted(model,
newdata = nd) %>%
data.frame() %>%
bind_cols(nd)
# Create new data for an average dog with a comorbidity
nd2 <- tibble(sFat = seq(from = -3.2, to = 3.2, length.out = 30), visit = "v0", breed = "mix", sex = "Female", comorbidity = "yes", sAge = 0)
# now use `fitted()` to get the model-implied trajectories
pr2 <- fitted(model,
newdata = nd2) %>%
data.frame() %>%
bind_cols(nd2)
# plot predictions for no comorbidity
ggplot(data = pr1, aes(x = sFat)) +
geom_smooth(aes(y = Estimate, ymin = Q2.5, ymax = Q97.5),
stat = "identity",
fill = "#F8766D", color = "#F8766D", alpha = 1/5, linewidth = 1.25) +
# plot predictions for comorbidity
geom_smooth(data = pr2, aes(y = Estimate, ymin = Q2.5, ymax = Q97.5),
stat = "identity",
fill = "#00BFC4", color = "#00BFC4", alpha = 1/5, linewidth = 1.25) +
# plot actual data
geom_point(data = df2,
aes(y = slgCort, colour = comorbidity),
size = 3) +
theme_bw() +
labs(title = "Predicted effect of body fat on log hair cortisol",
x = "Body fat (standardised)",
y = "Log hair cortisol (standardised)") +
theme(axis.title.y = element_text(size=14, face="bold"),
axis.title.x = element_text(size=14, face="bold"),
title = element_text(size=14, face="bold"),
plot.title = element_text(hjust = 0.5),
legend.position = "inside", legend.position.inside = c(0.87, 0.13)) +
coord_cartesian(ylim = c(-3, 3)) +
scale_x_continuous(breaks=seq(-3, 3 , 1))
##. b. Using ‘predict’ which is a wrapper for posterior_pred (returns error from model and residual error)
# Create new data for an average dog without a comorbidity
nd <- tibble(sFat = seq(from = -3.2, to = 3.2, length.out = 30), visit = "v0", breed = "mix", sex = "Female", comorbidity = "no", sAge = 0)
# now use `fitted()` to get the model-implied trajectories
pr1 <- predict(model,
newdata = nd) %>%
data.frame() %>%
bind_cols(nd)
# Create new data for an average dog with a comorbidity
nd2 <- tibble(sFat = seq(from = -3.2, to = 3.2, length.out = 30), visit = "v0", breed = "mix", sex = "Female", comorbidity = "yes", sAge = 0)
# now use `fitted()` to get the model-implied trajectories
pr2 <- predict(model,
newdata = nd2) %>%
data.frame() %>%
bind_cols(nd2)
# plot predictions for no comorbidity
ggplot(data = pr1, aes(x = sFat)) +
geom_smooth(aes(y = Estimate, ymin = Q2.5, ymax = Q97.5),
stat = "identity",
fill = "#F8766D", color = "#F8766D", alpha = 1/5, linewidth = 1.25) +
# plot predictions for comorbidity
geom_smooth(data = pr2, aes(y = Estimate, ymin = Q2.5, ymax = Q97.5),
stat = "identity",
fill = "#00BFC4", color = "#00BFC4", alpha = 1/5, linewidth = 1.25) +
# plot actual data
geom_point(data = df2,
aes(y = slgCort, colour = comorbidity),
size = 3) +
theme_bw() +
labs(title = "Predicted effect of body fat on log hair cortisol",
x = "Body fat (standardised)",
y = "Log hair cortisol (standardised)") +
theme(axis.title.y = element_text(size=14, face="bold"),
axis.title.x = element_text(size=14, face="bold"),
title = element_text(size=14, face="bold"),
plot.title = element_text(hjust = 0.5),
legend.position = "inside", legend.position.inside = c(0.87, 0.13)) +
coord_cartesian(ylim = c(-3, 4)) +
scale_x_continuous(breaks=seq(-3, 3 , 1))
draws <- as.matrix(model)
mean(draws[,3] >0)
## [1] 0.9865417
HPDI(draws[,3], prob=0.97)
## |0.97 0.97|
## 0.01884829 0.58929953
hypothesis(model, "sFat >0")
## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
## 1 (sFat) > 0 0.33 0.13 0.1 0.52 73.3 0.99 *
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
hypothesis(model, "sFat <0")
## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
## 1 (sFat) < 0 0.33 0.13 0.1 0.52 0.01 0.01
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
set.seed(666)
nd1 <- tibble(sFat = seq(-3, 3, by = 1), sAge = 0, visit = 0,
sex = "Female", comorbidity = "no",
breed = "mix",
case_number = c("dog1", "dog2", "dog3","dog4", "dog5", "dog6", "dog7"))
p1 <-
predict(model,
resp = "slgCort",
allow_new_levels = TRUE,
newdata = nd1) %>%
data.frame() %>%
bind_cols(nd1) %>%
ggplot(aes(x = sFat, y = Estimate, ymin = Q2.5, ymax = Q97.5)) +
geom_linerange(aes(ymin = Estimate - Est.Error, ymax = Estimate + Est.Error),
linewidth = 1, color = "#F8766D", alpha = 3/5) +
geom_point(size = 5, color = "#F8766D") +
theme_bw() +
labs(title = "Predicted effect of body fat on hair cortisol",
x = "Body fat (standardised)",
y = "Log hair cortisol (standardised)") +
theme(axis.title.y = element_text(size=12, face="bold"),
axis.title.x = element_text(size=12, face="bold"),
title = element_text(size=12, face="bold"),
plot.title = element_text(hjust = 0.5)) +
coord_cartesian(ylim = c(-2.5, 2.5)) +
scale_x_continuous(breaks=seq(-3, 3 , 1))
plot(p1)
set.seed(666)
nd3 <- tibble(sFat = seq(-3, 3, by = 1), sAge = 0, visit = 0,
sex = "Female", comorbidity = "no",
breed = "mix",
case_number = c("dog1", "dog2", "dog3","dog4", "dog5", "dog6", "dog7"))
nd4 <- tibble(sFat = seq(-3, 3, by = 1), sAge = 0, visit = 0,
sex = "Female", comorbidity = "yes",
breed = "mix",
case_number = c("dog1", "dog2", "dog3","dog4", "dog5", "dog6", "dog7"))
# predict outcome for comorbidity = no
pr3 <-
predict(model,
resp = "slgCort",
allow_new_levels = TRUE,
newdata = nd3) %>%
data.frame() %>%
bind_cols(nd3)
# predict outcome for comorbidity = yes
pr4 <-
predict(model,
resp = "slgCort",
allow_new_levels = TRUE,
newdata = nd4) %>%
data.frame() %>%
bind_cols(nd4)
# bind data for graphing
pr34 <-rbind(pr3, pr4)
pr34 <- as.data.frame(pr34)
p3 <- ggplot(pr34, aes(x = sFat, y = Estimate, , colour = comorbidity)) +
geom_pointrange(aes(ymin = Estimate - Est.Error, ymax = Estimate + Est.Error),
linewidth = 1, size = 0.5, position = position_dodge2(width = 0.5, preserve = "single")) +
theme_bw() +
labs(title = "Predicted effect of body fat & comorbidity on hair cortisol") +
labs(y = paste0("Predicted log hair cortisol (std)")) +
labs(x = paste0("Body fat (std)")) +
theme(axis.title.y = element_text(size=12, face="bold"),
axis.title.x = element_text(size=12, face="bold"),
title = element_text(size=12, face="bold"),
plot.title = element_text(hjust = 0.5)) +
coord_cartesian(ylim = c(-2.5, 2.5), xlim = c(-3.5, 3.5)) +
scale_x_continuous(limits = c(-3.5, 3.5), n.breaks = 7)
plot(p3)
# create new dataframe which contains results of the first dog
new_data <- rbind(df2[1,], df2[1,])
# Now change one category to be different
new_data$sFat <- c(-1, 1)
# Visualise df to make sure it has worked
new_data
## number group visit season breed_group coat_colour sex age comorbidity
## 1 c1 stopped v0 winter ret dark Male 43 yes
## 2 c1 stopped v0 winter ret dark Male 43 yes
## fat_percent cortisol lgCort breed sFat sAge slgCort
## 1 52.21393 4.92422 1.594166 ret -1 -1.485464 0.3415375
## 2 52.21393 4.92422 1.594166 ret 1 -1.485464 0.3415375
# Now get predictions from the draws of the model
pred_new <- posterior_predict(model, newdata = new_data)
# Compare difference in means between the two categories
difference <- pred_new[,2] - pred_new[,1]
# Examine mean of difference
mean(difference)
## [1] 0.6415713
# View histogram of this
dens(difference)
# Create HPDI
HPDI(difference, 0.97)
## |0.97 0.97|
## -3.096552 4.522566
# Create data frame for the plot
pred_new <- data.frame(pred_new)
pred_new$draws <- seq(1:nrow(pred_new))
colnames(pred_new) <- c("sFat-1", "sFat+1", "draws")
dens_l <- melt(pred_new, id.vars = "draws", variable_name = "standardised_fat")
# Plot distribution of difference
ggplot(dens_l, aes(value, fill = standardised_fat)) +
geom_density(alpha = 1/2) +
xlim(-5, 10) +
theme_bw() +
labs(title = "Predicted log hair cortisol by body fat mass") +
labs(y = paste0("Probability density")) +
labs(x = paste0("Predicted Log Hair Cortisol (standardised)")) +
theme(axis.title.y = element_text(size=12, face="bold"),
axis.title.x = element_text(size=12, face="bold"),
title = element_text(size=12, face="bold"),
plot.title = element_text(hjust = 0.5),
legend.position = "inside", legend.position.inside = c(0.87, 0.87))
# create new dataframe which contains results of the first dog
new_data2 <- rbind(df2[1,], df2[1,])
# Now change one category to be different
new_data2$sFat <- c(-2, 2)
# Visualise df to make sure it has worked
new_data2
## number group visit season breed_group coat_colour sex age comorbidity
## 1 c1 stopped v0 winter ret dark Male 43 yes
## 2 c1 stopped v0 winter ret dark Male 43 yes
## fat_percent cortisol lgCort breed sFat sAge slgCort
## 1 52.21393 4.92422 1.594166 ret -2 -1.485464 0.3415375
## 2 52.21393 4.92422 1.594166 ret 2 -1.485464 0.3415375
# Now get predictions from the draws of the model
pred_new2 <- posterior_predict(model, newdata = new_data2)
# Create data frame for the plot
pred_new2 <- data.frame(pred_new2)
pred_new2$draws <- seq(1:nrow(pred_new2))
colnames(pred_new2) <- c("sFat-2", "sFat+2", "draws")
dens_l2 <- melt(pred_new2, id.vars = "draws", variable_name = "standardised_fat")
# Plot distribution of difference
ggplot(dens_l2, aes(value, fill = standardised_fat)) +
geom_density(alpha = 1/2) +
xlim(-6, 10) +
theme_bw() +
labs(title = "Predicted log hair cortisol by body fat mass") +
labs(y = paste0("Probability density")) +
labs(x = paste0("Predicted Log Hair Cortisol (standardised)")) +
theme(axis.title.y = element_text(size=12, face="bold"),
axis.title.x = element_text(size=12, face="bold"),
title = element_text(size=12, face="bold"),
plot.title = element_text(hjust = 0.5),
legend.position = "inside", legend.position.inside = c(0.87, 0.87))
# create new dataframe which contains results of all dogs
new_data1 <- df2
# Now change one category to be different
new_data1$sFat <- -1
# create new dataframe which contains results of the first dog
new_data2 <- df2
# Now change one category to be different
new_data2$sFat <- 1
# Now get predictions from the draws of the models
pred_nd1 <- posterior_predict(model, newdata = new_data1)
pred_nd2 <- posterior_predict(model, newdata = new_data2)
pred_diff <- pred_nd2 - pred_nd1
pred_diff <- data.frame(pred_diff)
# Create mean of differences for each column (dog) of the dataframe
pred_diff_mean <- apply(pred_diff, 2, mean)
# View histogram of mean differences
dens(pred_diff_mean)
# mean difference
mean(pred_diff_mean)
## [1] 0.6508525
# Create HPDI
HPDI(pred_diff_mean, 0.97)
## |0.97 0.97|
## 0.6335775 0.6681280
# Set individual priors
prior_int <- set_prior("normal(0, 0.5)", class = "Intercept")
prior_b <- set_prior("normal(0, 1)", class = "b")
prior_b_sex <- set_prior("normal(0, 1)", class = "b", coef = "sexMale")
prior_b_co <- set_prior("normal(0.25, 1)", class = "b", coef = "comorbidityyes")
prior_b_f <- set_prior("normal(0, 0.5)", class = "b", coef = "sFat")
prior_b_sAge <- set_prior("normal(0, 0.5)", class = "b", coef = "sAge")
prior_sd <- set_prior("normal(0, 1)", class = "sd")
# Combine priors into list
priors2 <- c(prior_int, prior_b, prior_b_sex, prior_b_co, prior_b_f, prior_b_sAge, prior_sd)
Increased adapt_delta >0.8 (0.9999 here), as had divergent transitions
set.seed(666)
model2 <- brm(bf(slgCort ~ sFat + sAge + breed + sex + comorbidity + (1 | visit),
sigma ~ comorbidity),
family = gaussian(),
prior = priors2,
data = df2,
control=list(adapt_delta=0.9999, stepsize = 0.001, max_treedepth =15),
iter = 8000, warmup = 2000,
cores = 4,
save_pars = save_pars(all =TRUE),
sample_prior = TRUE)
## Compiling Stan program...
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 17.0.0 (clang-1700.0.13.5)’
## using SDK: ‘MacOSX15.5.sdk’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/opt/R/arm64/include -fPIC -falign-functions=64 -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
## 679 | #include <cmath>
## | ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1
## Start sampling
## Found more than one class "stanfit" in cache; using the first, from namespace 'rethinking'
## Also defined by 'rstan'
## Found more than one class "stanfit" in cache; using the first, from namespace 'rethinking'
## Also defined by 'rstan'
## Found more than one class "stanfit" in cache; using the first, from namespace 'rethinking'
## Also defined by 'rstan'
## Found more than one class "stanfit" in cache; using the first, from namespace 'rethinking'
## Also defined by 'rstan'
summary(model2)
## Family: gaussian
## Links: mu = identity; sigma = log
## Formula: slgCort ~ sFat + sAge + breed + sex + comorbidity + (1 | visit)
## sigma ~ comorbidity
## Data: df2 (Number of observations: 55)
## Draws: 4 chains, each with iter = 8000; warmup = 2000; thin = 1;
## total post-warmup draws = 24000
##
## Multilevel Hyperparameters:
## ~visit (Number of levels: 2)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.39 0.37 0.01 1.39 1.00 8922 9456
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept -0.80 0.45 -1.69 0.06 1.00 12410
## sigma_Intercept -0.90 0.33 -1.47 -0.18 1.00 12297
## sFat 0.25 0.13 -0.03 0.50 1.00 14898
## sAge 0.06 0.15 -0.24 0.36 1.00 21911
## breedckcs -0.06 0.50 -1.04 0.94 1.00 19989
## breedpug 0.11 0.48 -0.83 1.07 1.00 14806
## breedret -0.18 0.37 -0.90 0.55 1.00 14480
## breedother 0.39 0.39 -0.38 1.15 1.00 12968
## sexMale 0.22 0.28 -0.33 0.76 1.00 20566
## comorbidityyes 0.84 0.26 0.32 1.34 1.00 19495
## sigma_comorbidityyes 1.04 0.36 0.28 1.67 1.00 13377
## Tail_ESS
## Intercept 15550
## sigma_Intercept 13064
## sFat 14651
## sAge 17194
## breedckcs 16647
## breedpug 17360
## breedret 16487
## breedother 16384
## sexMale 16867
## comorbidityyes 16270
## sigma_comorbidityyes 14482
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
loo_model2 <- loo(model2, moment_match = TRUE)
loo_model2
##
## Computed from 24000 by 55 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -82.7 5.5
## p_loo 11.2 2.2
## looic 165.4 10.9
## ------
## MCSE of elpd_loo is 0.0.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 1.2]).
##
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
checks whether actual data is similar to simulated data.
color_scheme_set("blue")
pp_check(model2, ndraws = 100) +
theme_classic() +
labs(title = "Posterior predictions of body fat model (normal)",
y = "Density",
x = "Log hair cortisol (standardised") +
xlim(-5, 7) +
theme(axis.title.y = element_text(size=12, face="bold"),
axis.title.x = element_text(size=12, face="bold"),
title = element_text(size=12, face="bold"),
plot.title = element_text(hjust = 0.5),
axis.text.x = element_text(color = "grey50", size = 12),
axis.text.y = element_text(color = "grey8",size = 12), ,
legend.position = "inside", legend.position.inside = c(0.91, 0.83))
# Set individual priors
prior_int <- set_prior("normal(0, 0.5)", class = "Intercept")
prior_b <- set_prior("normal(0, 1)", class = "b")
prior_b_sex <- set_prior("normal(0, 1)", class = "b", coef = "sexMale")
prior_b_co <- set_prior("normal(0.25, 1)", class = "b", coef = "comorbidityyes")
prior_b_f <- set_prior("normal(0, 0.5)", class = "b", coef = "sFat")
prior_b_sAge <- set_prior("normal(0, 0.5)", class = "b", coef = "sAge")
prior_sd <- set_prior("normal(0, 1)", class = "sd")
# Combine priors into list
priors3 <- c(prior_int, prior_b, prior_b_sex, prior_b_co, prior_b_f, prior_b_sAge, prior_sd)
Increased adapt_delta >0.8 (0.9999 here), as had divergent transitions
set.seed(666)
model3 <- brm(bf(slgCort ~ sFat + sAge + breed + sex + comorbidity + (1 | visit),
sigma ~ comorbidity),
family = student(),
prior = priors3,
data = df2,
control=list(adapt_delta=0.9999, stepsize = 0.001, max_treedepth =15),
iter = 8000, warmup = 2000,
cores = 4,
save_pars = save_pars(all =TRUE),
sample_prior = TRUE)
## Compiling Stan program...
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 17.0.0 (clang-1700.0.13.5)’
## using SDK: ‘MacOSX15.5.sdk’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/opt/R/arm64/include -fPIC -falign-functions=64 -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
## 679 | #include <cmath>
## | ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1
## Start sampling
## Found more than one class "stanfit" in cache; using the first, from namespace 'rethinking'
## Also defined by 'rstan'
## Found more than one class "stanfit" in cache; using the first, from namespace 'rethinking'
## Also defined by 'rstan'
## Found more than one class "stanfit" in cache; using the first, from namespace 'rethinking'
## Also defined by 'rstan'
## Found more than one class "stanfit" in cache; using the first, from namespace 'rethinking'
## Also defined by 'rstan'
summary(model3)
## Family: student
## Links: mu = identity; sigma = log; nu = identity
## Formula: slgCort ~ sFat + sAge + breed + sex + comorbidity + (1 | visit)
## sigma ~ comorbidity
## Data: df2 (Number of observations: 55)
## Draws: 4 chains, each with iter = 8000; warmup = 2000; thin = 1;
## total post-warmup draws = 24000
##
## Multilevel Hyperparameters:
## ~visit (Number of levels: 2)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.39 0.38 0.01 1.42 1.00 10092 9944
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept -0.79 0.45 -1.67 0.09 1.00 13581
## sigma_Intercept -0.93 0.34 -1.53 -0.19 1.00 13282
## sFat 0.24 0.14 -0.05 0.49 1.00 14641
## sAge 0.04 0.16 -0.27 0.34 1.00 20236
## breedckcs -0.09 0.50 -1.07 0.89 1.00 20540
## breedpug 0.06 0.50 -0.90 1.08 1.00 15827
## breedret -0.18 0.37 -0.90 0.55 1.00 16563
## breedother 0.38 0.40 -0.40 1.16 1.00 13904
## sexMale 0.21 0.28 -0.34 0.75 1.00 21938
## comorbidityyes 0.82 0.27 0.28 1.33 1.00 17736
## sigma_comorbidityyes 1.01 0.37 0.23 1.69 1.00 14175
## Tail_ESS
## Intercept 15545
## sigma_Intercept 14262
## sFat 14787
## sAge 16248
## breedckcs 17568
## breedpug 16307
## breedret 17050
## breedother 16013
## sexMale 17538
## comorbidityyes 15625
## sigma_comorbidityyes 15339
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## nu 22.58 13.97 4.99 58.45 1.00 21794 14260
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
loo_model3 <- loo(model3, moment_match = TRUE)
loo_model3
##
## Computed from 24000 by 55 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -83.0 5.5
## p_loo 11.3 2.0
## looic 166.1 10.9
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 1.1]).
##
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
loo_compare(loo_model, loo_model2, loo_model3)
## elpd_diff se_diff
## model 0.0 0.0
## model2 -2.3 2.0
## model3 -2.7 2.0
checks whether actual data is similar to simulated data.
color_scheme_set("blue")
pp_check(model3, ndraws = 100) +
theme_classic() +
labs(title = "Posterior predictions of body fat model (Student's t)",
y = "Density",
x = "Log hair cortisol (standardised") +
xlim(-5, 7) +
theme(axis.title.y = element_text(size=12, face="bold"),
axis.title.x = element_text(size=12, face="bold"),
title = element_text(size=12, face="bold"),
plot.title = element_text(hjust = 0.5),
axis.text.x = element_text(color = "grey50", size = 12),
axis.text.y = element_text(color = "grey8",size = 12), ,
legend.position = "inside", legend.position.inside = c(0.91, 0.83))
# Set individual priors
prior_int <- set_prior("normal(0, 0.5)", class = "Intercept")
prior_sig <- set_prior("exponential(1)", class = "sigma")
prior_b <- set_prior("normal(0, 1)", class = "b")
prior_b_sex <- set_prior("normal(0, 1)", class = "b", coef = "sexMale")
prior_b_co <- set_prior("normal(0.25, 1)", class = "b", coef = "comorbidityyes")
prior_b_f <- set_prior("normal(0, 0.5)", class = "b", coef = "sFat")
prior_b_sAge <- set_prior("normal(0, 0.5)", class = "b", coef = "sAge")
prior_sd <- set_prior("normal(0, 1)", class = "sd")
prior_alpha <- set_prior("normal(4, 2)", class = "alpha")
# Combine priors into list
priors4 <- c(prior_int, prior_sig, prior_b, prior_b_sex, prior_b_co, prior_b_f, prior_b_sAge, prior_sd, prior_alpha)
Increased adapt_delta >0.8 (0.9999 here), as had divergent transitions
set.seed(666)
model4 <- brm(slgCort ~ sFat + sAge + breed + sex + comorbidity + (1 | visit),
family = skew_normal(),
prior = priors4,
data = df2,
control=list(adapt_delta=0.9999, stepsize = 0.001, max_treedepth =15),
iter = 8000, warmup = 2000,
cores = 4,
save_pars = save_pars(all =TRUE),
sample_prior = TRUE)
## Compiling Stan program...
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 17.0.0 (clang-1700.0.13.5)’
## using SDK: ‘MacOSX15.5.sdk’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/opt/R/arm64/include -fPIC -falign-functions=64 -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
## 679 | #include <cmath>
## | ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1
## Start sampling
## Found more than one class "stanfit" in cache; using the first, from namespace 'rethinking'
## Also defined by 'rstan'
## Found more than one class "stanfit" in cache; using the first, from namespace 'rethinking'
## Also defined by 'rstan'
## Found more than one class "stanfit" in cache; using the first, from namespace 'rethinking'
## Also defined by 'rstan'
## Found more than one class "stanfit" in cache; using the first, from namespace 'rethinking'
## Also defined by 'rstan'
summary(model4)
## Family: skew_normal
## Links: mu = identity; sigma = identity; alpha = identity
## Formula: slgCort ~ sFat + sAge + breed + sex + comorbidity + (1 | visit)
## Data: df2 (Number of observations: 55)
## Draws: 4 chains, each with iter = 8000; warmup = 2000; thin = 1;
## total post-warmup draws = 24000
##
## Multilevel Hyperparameters:
## ~visit (Number of levels: 2)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.40 0.38 0.01 1.40 1.00 8219 9090
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -0.27 0.50 -1.28 0.70 1.00 13870 14908
## sFat 0.18 0.17 -0.16 0.51 1.00 18969 17548
## sAge 0.01 0.15 -0.29 0.30 1.00 20284 16401
## breedckcs -0.12 0.41 -0.94 0.68 1.00 17551 16049
## breedpug 0.19 0.56 -0.92 1.26 1.00 17030 15958
## breedret -0.18 0.34 -0.83 0.49 1.00 15928 16747
## breedother 0.02 0.35 -0.65 0.73 1.00 14746 15376
## sexMale 0.15 0.25 -0.35 0.65 1.00 21009 16334
## comorbidityyes 0.40 0.36 -0.26 1.14 1.00 18046 15850
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.08 0.12 0.88 1.33 1.00 17460 16467
## alpha 3.93 1.63 1.04 7.47 1.00 14707 14861
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
loo_model4 <- loo(model4, moment_match = TRUE)
loo_model4
##
## Computed from 24000 by 55 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -85.1 5.7
## p_loo 9.7 2.3
## looic 170.1 11.3
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.6, 1.2]).
##
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
loo_compare(loo_model, loo_model4)
## elpd_diff se_diff
## model 0.0 0.0
## model4 -4.7 2.2
checks whether actual data is similar to simulated data.
color_scheme_set("blue")
pp_check(model4, ndraws = 100) +
theme_classic() +
labs(title = "Posterior predictions of body fat model (Student's t)",
y = "Density",
x = "Log hair cortisol (standardised") +
theme(axis.title.y = element_text(size=12, face="bold"),
axis.title.x = element_text(size=12, face="bold"),
title = element_text(size=12, face="bold"),
plot.title = element_text(hjust = 0.5),
axis.text.x = element_text(color = "grey50", size = 12),
axis.text.y = element_text(color = "grey8",size = 12), ,
legend.position = "inside", legend.position.inside = c(0.91, 0.83))